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We analyze the MeV /GeV emission from four bright Gamma- Ray Bursts (GRBs) observed by the 
Fermi- Large Area Telescope to produce robust, stringent constraints on a dependence of the speed of 
light in vacuo on the photon energy (vacuum dispersion) , a form of Lorentz invariance violation (LIV) 
allowed by some Quantum Gravity (QG) theories. First, we use three different and complementary 
techniques to constrain the total degree of dispersion observed in the data. Additionally, using a 
maximally conservative set of assumptions on possible source-intrinsic spectral-evolution effects, we 
constrain any vacuum dispersion solely attributed to LIV. We then derive limits on the “QG energy 
scale” (the energy scale that LIV-inducing QG effects become important, Eqg) and the coefficients of 
the Standard Model Extension. For the subluminal case (where high energy photons propagate more 
slowly than lower energy photons) and without taking into account any source-intrinsic dispersion, 
our most stringent limits (at 95% CL) are obtained from GRB 090510 and are Eqg,i > 7.6 times the 
Planck energy (Epf) and Eqg, 2 > 1.3 x 10 11 GeV for linear and quadratic leading order LIV-induced 
vacuum dispersion, respectively. These limits improve the latest constraints by Fermi and H.E.S.S. 
by a factor of ~ 2. Our results disfavor any class of models requiring Eqg,i < Ep\. 

PACS numbers: 11.30.Cp, 04.60.-m, 98.70.Rz 


I. INTRODUCTION 

While general relativity and Quantum Field Theory 
have each enjoyed impressive success so far, their for- 
mulations are currently inconsistent, hence motivating 
searches for unification schemes that can collectively be 
subsumed under the name of Quantum Gravity (QG) 
theories. These theories generally predict the existence 
of a natural scale at which the physics of space-time, 
as predicted by relativity theory, is expected to break 
down, hence requiring modifications or the creation of a 
new paradigm to avoid singularity problems. This scale, 
referred to as the “Quantum Gravity energy scale” Eqq, 
is in general expected to be of the order of the Planck 
scale [1], Ep\ = y 7 ( hc 5 )/G ~ 1.22 x 10 19 GeV, or in some 
cases lower (e.g., for some QG scenarios such as loop 
quantum gravity). 

Since relativity precludes an invariant length, the in- 
troduction of such a constant scale violates Lorentz In- 
variance (LI). Thus, tests of LI are strongly motivated 
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by the search for a theory of QG. Additional motivations 
for testing LI are the need to cut off high-energy (ultra- 
violet) divergences in quantum field theory and the need 
for a consistent theory of black holes [2, 3]. 

The idea that LI may be only approximate has been 
explored within the context of a wide variety of sug- 
gested Planck-scale physics scenarios. These include the 
concepts of deformed relativity, loop quantum gravity, 
non-commutative geometry, spin foam models, and some 
string theory (M theory) models (for reviews see, e.g., 
Refs. [4-6]). These theoretical explorations and their pos- 
sible consequences, such as observable modifications in 
the energy-momentum dispersion relations for free parti- 
cles and photons, have been discussed under the general 
heading of “Planck scale phenomenology”. 

There is also the motivation of testing LI in order to 
extent or limit its domain of applicability to the highest 
observable energies. Since the group of pure LI transfor- 
mations is unbounded at high energies, one should look 
for its breakdown at high energy scales, possibly through 
effects of Planck scale physics but perhaps through the 
effects of other unknown phenomena. To accomplish such 
a program, tests of the kinematics of LI violation (LIV) 
within the context of physical interaction dynamics such 
as quantum electrodynamics or standard model physics 
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(e.g., [7-11]) have been proposed. Fruitful frameworks for 
this kind of analysis, useful for testing the effects of LIV 
at energies well below the Planck scale, are the Taylor 
series expansion originally proposed in the seminal paper 
by Amelino-Camelia et al. [12] and the more compre- 
hensive Standard Model Extension (SME) parametriza- 
tion of Kostelecky and collaborators [13, 14]. These 
phenomenological parameterizations can be viewed as 
low-energy effective field theories, holding at energies 
E < Ep\ and providing an effective framework to search 
for LIV at energies far below the Planck scale. 

One manifestation of LIV is the existence of an energy- 
dependent “maximum attainable velocity” of a particle 
and its effect on the thresholds for various particle inter- 
actions, particle decays, and neutrino oscillations [7]. As- 
suming that the mass of the photon is zero, its maximum 
attainable velocity can be determined by measuring its 
velocity at the highest possible observable energy. This 
energy is, of necessity, in the gamma-ray range. Since 
we know that LI is accurate at accelerator energies, and 
even at cosmic-ray energies [15], any deviation of the ve- 
locity of a photon from its low energy value, c, must be 
very small at these energies. Thus, a sensitive test of LI 
requires high-energy photons (i.e., gamma rays) and en- 
tails searching for dependencies of the speed of light in 
vacuo on the photon energy. The method used in this 
study to search for such an energy dependence consists 
in comparing the time of flight between photons of differ- 
ent energies emitted together by a distant astrophysical 
source. As will be shown in the next section, the mag- 
nitude of a LIV-induced difference on the time of flight 
is predicted to be an increasing function of the photon 
energy and the distance of source. Thus, because of the 
high-energy extend of their emission (up to tens of GeV) , 
their large distances (redshifts up to a value of ~8), and 
their rapid (down to ms scale) variabilities, Gamma-Ray 
Bursts (GRBs) are very effective probes for searching for 
such LIV-induced spectral dispersions [12], 

There have been several searches for LIV applying 
a variety of analysis techniques on GRB observations. 
Some of the pr e-Fermi studies are those by Lamon et 
al. [16] using INTEGRAL GRBs; by Bolmont et al. [17] 
using HETE-2 GRBs; by Ellis et al. [18] using HETE, 
BATSE, and Swift GRBs; and by Rodrfguez-Martfnez 
et al. using Swift and Konus-Wind observations of 
GRB 051221 A [19]. The most stringent constraints, how- 
ever, have been placed using Fermi observations, mainly 
thanks to the unprecedented sensitivity for detecting the 
prompt MeV/GeV GRB emission by the Fermi Large 
Area Telescope (LAT) [20, 21], These constraints include 
those by the Fermi LAT and Gamma-Ray Burst Moni- 
tor (GBM) collaborations using GRBs 080916C [22] and 
090510 [23], and by Shao et al. [24] and Nemiroff et al. [25] 
using multiple Fermi GRBs. In addition to these GRB- 
based studies, there have been some results using TeV ob- 
servations of bright Active Galactic Nuclei flares, includ- 
ing the MAGIC analysis of the flares of Mrk 501 [26, 27] 
and the H.E.S.S. analysis on the exceptional flare of 


PKS 2155-304 [28, 29]. It should be noted that the stud- 
ies above did not assume any dependence of LIV on the 
polarization of the photons, manifesting as birefringence. 
In the case that such a dependence exists, constraints on 
LIV effects can be produced [10, 30, 31] that are 13 or- 
ders of magnitude stronger than the dispersion-only con- 
straints placed with time-of-flight considerations (as in 
this work). It should be added that there is a class of 
theories that allow for photon dispersion without bire- 
fringence that can be directly constrained by our results 
(e.g., [32]). 

The aim of this study is to produce a robust and com- 
petitive constraint on the dependence of the velocity of 
light in vacuo on its energy. Our analysis is performed 
on a selection of Fermi - LAT [21] GRBs with measured 
redshifts and bright GeV emission. We first apply three 
different analysis techniques to constrain the total degree 
of spectral dispersion observed in the data. Then, using 
a set of maximally conservative assumptions on the pos- 
sible source-intrinsic spectral evolution (which can mas- 
querade as LIV dispersion), we produce constraints on 
the degree of LIV-induced spectral dispersion. The lat- 
ter constraints are weaker than those on the total degree 
of dispersion, yet considerably more robust with respect 
to the presence of a source-intrinsic effects. Finally, we 
convert our constraints to limits on LlV-model-specific 
quantities, such as Eqq and the coefficients of the SME. 

The first method used to constrain the degree of dis- 
persion in the data, named “PairView” (PV) , is created as 
part of this study, and performs a statistical analysis on 
all the pairs of photons in the data to find a common 
spectral lag. The second method, named “Sharpness- 
Maximization Method” (SMM), is a modification of ex- 
isting techniques (e.g., DisCan [33]) and is based on the 
fact that any spectral dispersion will smear the struc- 
ture of the light curve, reducing its sharpness. SMM’s 
best estimate is equal to the negative of a trial degree of 
dispersion that, when applied to the actual light curve, 
restores its assumed-as-initially-maximal sharpness. Fi- 
nally, the third method employs an unbinned maximum 
likelihood (ML) analysis to compare the data as observed 
at energies low enough for the LIV delays to be negligible 
to the data at higher energies. The three methods were 
tested using an extensive set of simulations and cross- 
checks, described in several appendices. 

Our constraints apply only to classes of LIV models 
that possess the following properties. First, the magni- 
tude of the LIV-induced time delay depends either lin- 
early or quadratically on the photon energy. Second, 
this dependence is deterministic, i.e., the degree of LIV- 
induced increase or decrease in the photon propagation 
speed does not have a stochastic (or “fuzzy”) nature as 
postulated by some of the LI models (see Ref. [34] and 
references therein). Finally, the sign of the effect does 
not depend on the photon polarization - the velocities 
of all photons of the same energy are either increased or 
decreased due to LIV by the same exact amount. 

In Sec. II we describe the LIV formalism and nota- 
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tion used in the paper, in Sec. Ill we describe the data 
sets used for the analysis, in Sec. IV we describe the 
three analysis methods and the procedure we used to take 
into account possible intrinsic spectral-evolution effects, 
in Sec. V we present the results, in Sec. VI we report 
and discuss their associated systematic uncertainties and 
caveats, and finally in Sec. VII we compare our results to 
previous measurements and discuss their physical impli- 
cations. We present our Monte Carlo simulations used for 
verifying PV and SMM in Appendix A, the calibrations 
and verification tests of the ML method in Appendix B, 
a direct one-to-one comparison of the results of the three 
methods after their application on a common set of simu- 
lated data in Appendix C, and cross-checks of the results 
in Appendix D. 


II. FORMALISM 

In QG scenarios, the LIV-induced modifications to the 
photon dispersion relation can be described using a series 
expansion in the form 


E 2 ~ p 2 c 2 x 



where c is the constant speed of light (at the limit of 
zero photon energy), s± is the “sign of LIV”, a theory- 
dependent factor equal to +1 (—1) for a decrease (in- 
crease) in photon speed with an increasing photon energy 
(also referred to as the “subluminal” and “super luminal” 
cases). For E <C Eqq , the lowest order term in the se- 
ries not suppressed by theory is expected to dominate 
the sum. In case the n = 1 term is suppressed, some- 
thing that can happen if a symmetry law is involved, the 
next term n = 2 will dominate. We note that in effective 
field theory n = d — 4, where d is the mass dimension of 
the dominant operator. Therefore, the n = 1 term arises 
from a dimension 5 operator [35] . It has been shown that 
odd mass-dimension terms violate CVT [13, 36]. Thus, 
if CVT is preserved, then the n = 2 term is expected to 
dominate. In this study, we only consider the n = 1 and 
n = 2 cases, since the Fermi results are not sensitive to 
higher order terms. 

Using Eq. 1 and keeping only the lowest-order domi- 
nant term, it can be found that the photon propagation 
speed itph, given by its group velocity, is 


Uph(E) 




n + 1 
2 



(2) 


where c = bin u p h(E). Because of the energy dependence 

of Uph(E), two photons of different energies Eh > E\ 
emitted by a distant source at the same time and from 
the same location will arrive on Earth with a time delay 
At which depends on their energies. We define the “LIV 


parameter” r„ as the ratio of this delay over E^-E™ [37]: 
At (1 + n) 1 
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is a comoving distance that also depends on the order of 
LIV (n), z is redshift, Hq is the Hubble constant, and 
Ha and SIm are the cosmological constant and the total 
matter density (parameters of the ACDM model). 

In the SME framework [14], the slight modifications in- 
duced by LIV effects are also described by a series expan- 
sion with respect to powers of the photon energy. In this 
framework, LIV can also be dependent on the direction 
of the source. Including only the single term assumed to 
dominate the sum, r n is given by: 


Tn - H 0*Wn)c["j2 X Kn, (5) 

\ira / 

where n is the direction of the source, oXjrn(iT.) are spin- 
weighted spherical harmonics, and are coefficients 

of the framework that describe the strength of LIV. In the 
SME case of a direction-dependent LIV, we constrain the 
sum (enclosed in parentheses) as a whole. For the alter- 
native possibility of direction independence, all t he terms 
in the sum become zero except o^oo = Voo = -\/l/(47r). 
In that case, we constrain a single coefficient. 

The coordinates of n are in a Sun-centered celestial 
equatorial frame described in Section V of Ref. [14], di- 
rectly related to the equatorial coordinates of the source 
such that its first coordinate is equal to 90° — Declination 
and the second being equal to the Right Ascension. Fi- 
nally, the coefficients can be either positive or neg- 

ative depending on whether LIV-induced dispersion cor- 
responds to a decrease or increase in photon speed with 
an increasing energy respectively (i.e., the sign of the 
SME coefficients plays the role of the s± factor of the 
series-expansion framework) . 

In the important case of a d = 5 modification of the 
free photon Lagrangian in effective field theory, Myers 
and Pospelov have shown that the only d = 5 (n = 1) 
operator that preserves both gauge invariance and rota- 
tional symmetry implies vacuum birefringence [35]. In 
such a case, and as was mentioned in the Introduction, 
significantly stronger constraints can be placed using the 
existence of birefringence than with just dispersion (as in 
this work). For this reason, in this paper and when work- 
ing within the given assumptions of the SME framework, 
we proceed assuming that the d = 5 terms are either zero 
or dominated by the higher-order terms, and proceed to 
constrain the d = 6 terms, which are not expected to 
come with birefringence. 
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Our aim is to constrain the LIV-related parameters in- 
volved in the above two parametrizations: the quantum 
gravity energy Eqq (for n = {1,2} and s± = ±1) and 
the coefficients of the SME framework (for d = 6 for 
both the direction dependent and independent cases) . 
To accomplish this, we first constrain the total degree of 
spectral dispersion in the data, r n , and then using the 
measured distance of the GRB and the cosmological con- 
stants, we calculate lower limits on Eqq through Eq. 3 
and confidence intervals for the SME coefficients using 
Eq. 5. We also produce an additional set of constraints 
after accounting for GRB-intrinsic spectral evolution ef- 
fects (which can masquerade as LIV). In that case, we 
first treat r n as being the sum of the GRB-intrinsic dis- 
persions Ti n t and the LIV-induced dispersion tliv, then 
we constrain tliv assuming a model for Tint, and finally 
constrain Eqq and the SME coefficients using the con- 
straints on tliv. 

We employ the cosmological parameters determined 
using WMAP 7-year data Hm = 0.272 and = 
0.728 [38], and a value of Ho = 73.8 ± 2.4 kms -1 Mpc -1 
as measured by the Hubble Space Telescope [39]. 

III. THE DATA 

We analyze the data from the four Fermi - LAT 
GRBs having bright GeV prompt emission and 
measured redshifts, namely GRBs 080916C, 090510, 
090902B, and 090926 A. We analyze events passing the 
P7_TRANSIENT_V6 selection, optimized to provide in- 
creased statistics for signal- limited analyses [40]. Its 
main difference from the earlier P6_V3_TRANSIENT se- 
lection used to produce previous Fermi constraints on 
LIV [22, 23] consists in improvements in the classification 
algorithms, which brought an increase in the instrument’s 
acceptance mostly below ^300 MeV 1 . 

We reject events with reconstructed energies less than 
30 MeV because of their limited energy and angular re- 
construction accuracy. We do not apply a maximum- 
energy cut. In the case of GRB 080916C, however, we 
removed an 106 GeV event detected during the prompt 
emission, since detailed analyses by the LAT collabora- 
tion 2 showed that it was actually a cosmic-ray event mis- 
classified as a photon. 

We keep events reconstructed within a circular region 
of interest (ROI) centered on the GRB direction and of a 
radius large enough to accept 95% of the GRB events ac- 
cording to the LAT instrument response functions, i.e. , a 
radius equal to the 95% containment radius of the LAT 
point spread function (PSF). Because the LAT PSF is 
a function of the true photon energy and off-axis angle 
(the angle between the photon true incoming direction 


1 For a detailed list of differences see http://fermi.gsfc.nasa.gov/ 
ssc/data/analysis/documentation/Pass7_usage.html 

2 not published 


TABLE I. Distances of Analyzed GRBs 


GRB 

Redshift 

ft 1 

ft 2 

080916C 

4.35 ± 0.15 [41] 

4.44 

13.50 

090510 

0.903 ± 0.003 [42] 

1.03 

1.50 

090902B 

1.822 [43] a 

2.07 

3.96 

090926A 

2.1071 ±0.0001 [44] 

2.37 

4.85 


This GRB had a spectroscopically-measured redshift, which 
implies an error at the 10 — 3 level. 


and the LAT boresight), the PSF containment radius is 
calculated on a per-photon basis. In this calculation, we 
approximate the (unknown) true off-axis angles and en- 
ergies with their reconstructed values, something that in- 
duces a slight error at low energies. Below ^100 MeV, the 
LAT angular reconstruction accuracy deteriorates and 
the 95% containment radius becomes very large. To limit 
the inclusion of background events due to a very large 
ROI radius and also reject some of the least accurately 
reconstructed events, we limit the ROI radius to be less 
than 12°. The GRB direction used for the ROI’s center 
is obtained by follow-up ground-based observations (see 
citations in Tab. I) and can be practically assumed to 
coincide with the true direction of the GRB. 

The above data set are further split and cut depending 
on the requirements of each of the three analysis methods 
(as described below). Figure 1 shows the light curves and 
the event time versus energy scatter plots of the GRBs 
in our sample, and Tab. I shows the GRB redshifts and 
K\ and K 2 distances (defined in Eq. 4). 

Time Interval Selection 

The analyzed time intervals are chosen to correspond 
to the period with the highest temporal variability, fo- 
cusing on the brightest pulse of each GRB. This choice 
is dictated by the fact that GRB emission typically ex- 
hibits spectral variability, which can potentially manifest 
as a LIV-dispersion effect (see discussion in Sec. VI for 
details on GRB spectral variability). By focusing on a 
narrow snapshot of the burst’s emission, we aim to ob- 
tain constraints that are affected as little as possible by 
such GRB-intrinsic effects. Starting from this require- 
ment, we select the time intervals to analyze, hereafter 
referred to as the “default” time intervals, using a pro- 
cedure we devised a priori and applied identically on all 
four GRBs. 

We start by characterizing the brightest pulse in each 
GRB by fitting its time profile with the flexible model 
used by Norris et al. [45] to successfully fit more than 
400 pulses of bright BATSE bursts: 

Aexp[ (|t tmax | / @r) ] 1 ^ ^max ( /^\ 

Aexp[-(|t - t ma , x \/a d ) v ] t > f max , 

where f max is the time of the pulse’s maximum intensity 
A, v is a parameter that controls the shape of the pulse, 
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FIG. 1. Time and energy profiles of the detected events from the four GRBs in our sample. Each column shows an event energy 
versus event time scatter plot (top) and a light curve (bottom). The vertical lines denote the time intervals analyzed (solid line 
for n = 1 and dashed line for n = 2), the choice of which is described in Sec. IV. If a dashed line is not visible, it approximately 
coincides with the solid one. 


and ay and oy are the rise and decay time constants. 
For v = {1,2} the equation describes a two-sided expo- 
nential or Gaussian function respectively. We use the 
best fit parameters (as obtained from a maximum likeli- 
hood analysis) to define a “pulse interval” extending from 
the time instant that the pulse height rises to 5% of its 
amplitude to the time instant that it fells to 15% of its 
amplitude. We choose such an asymmetric cut because 
of the long falling-side tails of GRB pulses. 

We then expand this initial “pulse interval” until no 
photons that were generated outside of it (at the source) 
could have been detected inside of it (at the Earth) due 
to LIV dispersion, and also until no photons that were 
generated inside of it (at the source) could have been 
detected outside of it (at the Earth) due to LIV disper- 
sion. We use conservative values of Eqg,i = 0.5 x Epi 
and Eqq 2 = 1.5 x 10 10 GeV for the maximum degree 
of LIV dispersion considered in extending the time inter- 
val, values which correspond to roughly one half of the 


stringent and robust limits obtained by Fermi [23] and 
H.E.S.S. [28, 29]. The interval resulting from this expan- 
sion is the one chosen for the analysis (hereafter referred 
to as the “default” interval). The main reason for ex- 
tending the interval is to avoid constraining the possible 
emission time of the highest-energy photons in the initial 
“pulse interval” to a degree that would imply an artifi- 
cially small level of dispersion. 


The choice of time interval for GRB 090510 and n = 1 
is demonstrated in Fig. 2. The (default) time intervals 
for all GRBs are shown in Fig. 1 with the vertical solid 
(n = 1) and dashed lines (n = 2), and are also reported 
in Tab. II. 
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FIG. 2. Demonstration of the calculation of the default time 
interval for GRB 090510 and linear LIV. The data points show 
the event rate for energies greater than 30 MeV focusing on its 
main pulse, the thick curve shows the fit on the pulse using 
Eq. 6, the dashed vertical lines denote the “pulse interval”, 
and the solid vertical lines denote the extended time interval 
(the default interval) chosen for the analysis. 


IV. DATA ANALYSIS 

A. PairView and Sharpness-Maximization Methods 

Because the way we calculate confidence intervals is 
identical between PV and SMM, we first describe how 
the best estimate of the LIV parameter is calculated by 
each of these two methods, and then proceed to describe 
their common confidence-interval calculation procedure. 


Best-Estimate Calculation: PairView 

The PV method calculates the spectral lags l t j be- 
tween all pairs of photons in a data set and uses the dis- 
tribution of their values to estimate the LIV parameter. 
Specifically, for a data set consisting of N photons with 
detection times and energies Ld .jv, the method 

starts by calculating the N x (AT — 1)/2 photon-pair spec- 
tral lags l t ,j for each i > j: 

k i = U ~ tj (7) 

(where n is the order of LIV), and creates a distribution 
of their values. 

Let us examine how the distribution of lij values de- 
pends on the properties of the data and LIV disper- 
sion. For a light curve comprising at emission a single 
5- function pulse and for a dispersion r„, the lij distri- 
bution will consist of a single (5-function peak at a value 
of exactly r„. For a light curve comprising (at emission) 
a finite-width pulse, the now non-zero time differences 
between the emission times of the events behave as noise 


inducing a non-zero width to the distribution of lij . Sim- 
ilarly to the previous ideal case however, the Uj distribu- 
tion will be peaked at approximately r n . For a realistic 
light curve consisting of one or more peaks superimposed 
on a smoothly varying emission, the distribution of lij 
will be composed of a signal peak centered at ~ r n (con- 
sisting of lij values created primarily by events i, j emit- 
ted temporally close and with not too similar energies) 
and a smoother underlying wide background (consisting 
of the rest of the kj values). 

Following the above picture, the estimator f n of r„ is 
taken as the location of the most prominent peak in the 
lij distribution. This peak becomes taller and narrower, 
thus more easily detectable, as the variability time scale 
decreases and as the width of the energy range increases. 

Searching for the peak using a histogram of the lij 
values would require us to first bin the data, a procedure 
that would include choosing a bin width fine enough to 
allow for identifying the peak with good sensitivity but 
also wide enough to allow for good statistical accuracy 
in the bin contents. We decided not to use a histogram 
to avoid the subjective choice of bin width. Instead, we 
use a kernel density estimate (KDE), as it provides a way 
to perform peak finding on unbinned data, and as it is 
readily implemented in easy to use tools with the ROOT 
TKDE method 3 . We use a Gaussian kernel for the KDE 
and a bandwidth chosen such as to minimize the Mean 
Integrated Squared Error calculated between the KDE 
and a very finely binned histogram of the photon-pair 
lags. 

Best-Estimate Calculation: Sharpness- Maximization Method 

SMM is based on the fact that the application of any 
form of spectral dispersion to the data (e.g., by LIV) will 
smear the light curve decreasing its sharpness. Based 
on this, SMM tries to identify the degree of dispersion 
that when removed from the data (i.e. , when the negative 
value of it is applied to the data) maximizes its sharp- 
ness. This approach is similar to the “Dispersion Cancel- 
lation” (DisCan) technique [33] , the “Minimal Dispersion” 
method [46], and the “Energy Cost Function” method 
[26, 46]. The most important difference between these 
approaches is the way the sharpness of the light curve is 
measured. 

We start the application of SMM by analyzing a data 
set consisting of photons with detection times ti and en- 
ergies Ei to produce a collection of “inversely smeared” 
data sets, each corresponding to a trial LIV parameter 
T n , by subtracting E ] 1 x r„ from the detection times ti. 
For each of the resulting data sets, the modified pho- 
ton detection times are first sorted to create a new set 
and then the sharpness of its light curve is measured 


3 http://root.cern.ch/root/html/TKDE.html 
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using t\ and Ei. After this procedure has been applied 
on a range of trial LIV parameters, we find the inversely 
smeared data set with the sharpest light curve, and select 
the trial r n value used to produce it as the best estimate 

of T n . 

In their analysis of the data from a flare of the blazar 
Mrk 501, the MAGIC collaboration [26] quantified the 
sharpness of the light curve using an “Energy Cost Func- 
tion”, which was essentially the sum of the photon ener- 
gies detected in some predefined time interval chosen to 
correspond to the most active part of the flare. Scargle et 
al. [33] explored a range of different cost functions to mea- 
sure the sharpness of the light curve, including Shannon, 
Renyi, and Fisher information, variance, total variation, 
and self-entropy, finding that the Shannon information 
is the most sensitive. In this study, we use a function S 
that is similar to the Shannon information and is defined 
as: 

N-p / \ 

5 ( r ") = 51 lo s f ’ ( 8 ) 

i= 1 \ i+P V 

where p is a configurable parameter of the method. 

Different values of p will tune the algorithm to evaluate 
the sharpness of the light curve focusing on intervals con- 
sisting of different numbers of events (i.e., of p events) or 
equivalently focusing on different time scales. As a result, 
the choice of p affects the performance of the algorithm in 
two ways. For a small value of p (up to ^3), some of the 
durations in the denominator of S(r n ) can become rela- 
tively very small, making some of the 1 / (t' i+p — t'fi) terms 
very large. In this case, S(r n ) can fluctuate significantly 
as a function of the trial lag, decreasing the accuracy with 
which the best LIV parameter can be measured. For too 
large values of p the algorithm essentially tries to mini- 
mize the total duration of the analyzed data, focusing on 
time scales larger than the variability time scale, ending 
up with a diminished sensitivity (in practice the peak of 
S becomes flatter). These effects are demonstrated in 
Fig. 3. 

To choose the value of p we first generate a large num- 
ber of simulated data sets inspired by the GRB under 
study, we then apply the method using a series of dif- 
ferent p values, and finally we choose the p value that 
produces the most constraining median upper limit on r„ 
(for s± = +1). These simulated data sets are constructed 
similarly to the procedure described in Appendix A us- 
ing a light-curve template produced by a KDE of the 
actual light curve, with the same statistics as the data, 
a spectrum similar to that in the data, and without any 
spectral dispersion applied. 

Finally, it should be noted that the method’s descrip- 
tion above was for the case of zero source-intrinsic spec- 
tral evolution effects since the light curve of the GRB 
mission at the source was treated as being maximally 
sharp. This picture is equivalent to assuming that there 
is an initial (imaginary) maximally-sharp signal that is 
first distorted by GRB-intrinsic effects and then by LIV. 


In that case, the constraints provided by SMM will be on 
the aggregate effect. 



FIG. 3. Curves of SMM’s sharpness measure <S(r„) versus 
the trial value of the LIV parameter n, each produced using 
a different p value. These curves were generated using the 
GRB 090510-inspired data set described in Appendix A after 
the application of a dispersion equal to +0.04 s/GeV (value 
denoted with the vertical line) . The circles denote the maxima 
of the curves, the positions of which are used to produce T\ . 
As can be seen, too small or too large values of p correspond 
to a reduced accuracy for measuring the position of the peak. 


Confidence-Interval Calculation 

The PV and SMM methods produce a confidence in- 
terval on the best LIV parameter by means of a random- 
ization analysis. 

We start by producing one hundred thousand random- 
ized data sets by shuffling the association between ener- 
gies and times of the detected events. Because the total 
number of events and the distributions of energies and 
times are identical between the actually detected and the 
randomized data sets, their statistical power (i.e., their 
ability to constrain the dispersion) is similar. However, 
because of the randomization, any dispersion potentially 
present in the actual data is lost. After the set of ran- 
domized data sets is constructed, the best LIV parameter 
is measured on each one of them and the measurements 
are used to create a (normalized to unity) distribution 

fr- 

We then define the measurement error on r n (for the 
general case of any r n ) as £ = f n — r n and the probability 
distribution function (PDF) of £ as Ps(e), where e is a 
random realization of £ . We assume that Ps has a negli- 
gible dependence on r n (at least for the range of values of 
r n expected to be present in the data) and approximate: 

P e (e)~P £ (e\Tn = 0). (9) 

The PDF Ps for the case of a zero r„, Pg(e|r n = 0), can 
be identified as the normalized distribution f r produced 


using our randomization simulations. Thus, 


P £ (e) ~ f r (e). (10) 

Since £ is a quantity with a known PDF and since it 
depends on the unknown parameter r n , it can be used as 
a pivotal quantity to construct a two-sided Confidence 
Interval (Cl) of Confidence Level (CL) for r„ as: 

CL = Pr(q(i _ CL )/2 <£ < q(i+CL)/2) (H) 

= P r (<I(l-CL)/2 < T n — T n < <7(i+CL)/2) 

= Pr{T n — q(l+CL)/2 < T n < T n — q(l-CL)/2) 

= Pr(LL <Tn< UL), 

where LL = r n - q(i+cL )/2 and UL = r n - q(i-CL )/2 are 
the lower and upper limits defining the Cl, and < 7 (i_cl )/2 
and q(i+cL )/2 are the (1 — CL)/2 and ( l+CL)/2 quantiles 
of f r . 

To produce a lower limit on Eqq for the subluminal 
or the superluminal case, we use eq. 3 substituting r n 
with its lower or upper limit, respectively, and solve for 
pQG- 


B. Likelihood Method 

The ML fit procedure used in this work has been de- 
veloped and applied by Martinez and Errando [27] to 
MAGIC data for the 2005 flare of Mkn 501 and by 
Abramowski et al. [29] to H.E.S.S. data for the gigantic 
flare of PKS 2155—304 in 2006. This section describes its 
key aspects, its underlying assumptions, and the details 
of its application to GRB data. 

The ML method consists in comparing the arrival time 
of each detected photon with a template light curve which 
is shifted in time by an amount depending linearly or 
quadratically on the event’s energy. For a fixed num- 
ber of independent events TVgt with energies and times 
{Ei, ti}i=i,jVf- lt observed in the energy and time intervals 
[E cut ,E max ] and [t i , ^ 2 ] , the unbinned likelihood function 
is: 


N iit 

£ = P(Ei,ti\T n ), 

i— 1 


( 12 ) 


where P is the PDF of observing one event at en- 
ergy E and time f, given r n . For an astrophysi- 
cal source observed by a gamma-ray telescope, it is 
P(Ei,U\T n ) = R(Ei,ti\T n )/N pied , where R is the ex- 
pected differential count rate at energy E and time t and 

pE max pt2 

Npred = / R(E,t\r n ) cLE dt is the total number 

J Ec ut J t\ 

of events predicted by the model. For a point-like source 
observed by the Fermi- LAT: 

/>00 

R(E,t\T n )= / F(E t ,t\T n ) A eS (E t ) T>(E t ,E) dE t , 

Jo 

( 13 ) 


where F(E t ,t\T n ) is the model for the photon flux which 
is incident on the LAT at the photon (true) energy E t 
and time t, whereas A e s(E t ) 4 and V(E t , E) are the LAT 
effective area and energy redistribution functions, respec- 
tively. As the energy resolution with the LAT is better 
than 15% above 100 MeV [40] we can neglect any energy 
mis-reconstruction effects. Assuming no spectral vari- 
ability and that the flux spectrum follows a power law 
with possible attenuation at the highest energies, then: 

F(E, t\r n ) = 0o E~ T e- E ' Ei f{t - T n E n ), (14) 

where F is the time-independent spectral index, Ef is the 
cutoff energy, and the function f(t) is the time profile 
of the emission that would be received by the LAT in 
case of a null LIV-induced lag T n E n . We explain further 
below how the function f[t) is derived from the data in 
practice. Finally, defining the observed spectral profile 
as A (E) = 0o E~ r e~ E ! E * A e ${E), we obtain: 

P{E, t\r n ) = A [E) fit - T n E n )/N pied (15) 


. Thus, the ML estimator r n of the LIV parameter r n 
satisfies: 


y- 9 log f(U - r n E?) 


N fit cWpred 
N pred Jt T : 


T n =Tn 


= 0 . 


(16) 

For the brightest LAT-detected GRBs, iVg t > 50 typi- 
cally (see Tab. II) thus a good estimate of the sensitivity 
offered by the estimator f n can be obtained by consid- 
ering the ideal case of the large sample limit. In this 
regime, f n is unbiased and efficient like any ML estima- 
tor. Namely, its variance reaches the Cramer-Rao bound, 
e.g., given by Eq. (9.34) page 217 of [47]: 



As the time profile can be measured up to very large times 
in case of large photon statistics, one can show that the 
standard deviation of f n is simply given by: 




(18) 


/ -+- OO 

fit) 2 /fit) dt (= 1 / cr 2 for a Gaus- 

-OO 

sian time profile of standard deviation cr), ( E mxn )h = 


E mxn A iE) dE / A h and A h = 


A iE) dE. 


4 The effective area also depends on the direction of the source in 
instrument coordinates, a typically continuously varying quan- 
tity. We can drop the time dependence by approximating 
A e ff (Et, t) with its averaged over the observation value A g r (E t). 
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The above expression for cr[f n ] is a good approximation 
(within a factor 2 to 3) of the actual standard deviation 
of T n , and it gives a useful estimate of the expected sen- 
sitivity. However, our final results are based on a proper 
derivation of confidence intervals as described further be- 
low in this section. 

The spectral profile A (E) is constant with time since 
T is assumed to be constant during the considered time 
interval (see further discussions on possible spectral evo- 
lution effects in Sec. VI). The spectral profile is also in- 
dependent of the LIV parameter, and is only used as 
a weighting function in the PDF normalization V pre( j. 
For these reasons, we approximate the spectral profile 
by a power-law function (with a fixed attenuation when 
needed) : 


A(E) oc FT 7 e~ E/Ei . (19) 

The spectral index 7 is obtained from the fit of the above 
function to the time- integrated spectrum S(E) observed 
by the LAT: 


rt 2 


S(E)= F{E,t\T n )A eS {E)dt. 


(20) 


In practice, we define a E cu t £ [100, 150] MeV and we 
fit the spectrum S(E) above E cut (see Fig. 7 for an ex- 
ample) in order to obtain a fairly good estimate of the 
spectral index, namely 7 ~ F within errors (see discus- 
sion in Sec. VI regarding this approximation). For the 
case of GRB 090926A, we use a power-law function that 
has an exponential break, in accordance with the findings 
of Ackermann et al. [48]. 

Knowledge of the time profile f(t) is crucial for the ML 
analysis. Typically, E cu t divides the LAT data set in two 
samples of roughly equal statistics. The ML analysis is 
performed using events with energies above E cnt , whereas 
the fit of the light curve C ( t ) observed by the LAT below 
E cu t is used to derive the time profile: 


C{t) = / A (E) f(t - r n E n ) dE ~ A; f[t - r n <£"),], 

“ -E'min 

( 21 ) 

r E cut 

A (E) dE and {E n )i = 


where E m [ n = 30 MeV, A ; = 

pE cut 




E n A (E) dE / A/. The Taylor expansion used in 


'E a 


Eq. (21) is justified as LIV-induced lags are effectively 
negligible for low-energy events, and it yields the time 
profile: 


for better statistics and because the calculations need an 
estimate of the GRB flux at times that are also external 
to the default time intervals. 

We then proceed with calculating the likeli- 
hood function C for a series of trial values of 
the LIV parameter 77, and plotting the curve of 
— 2Aln(£) = — 21n [£(r n )/£(f n )] as a function of r n . We 
first produce a best estimate of r„, t„, equal to the 
location of the minimum of the — 2Aln(£) curve. We 
also produce a Cl on r n for an approximately two-sided 
CL (90%, 99%) using the two values of T n around the 
global minimum at f n for which the curve reaches a 
values of 2.71 and 6.63, respectively 5 . Hereafter we 
refer to these CIs as being obtained “directly from the 
data”. In addition, we produce a set of “calibrated” CIs 
on t n using Monte Carlo simulations and as described 
in Appendix B. The calibrated CIs take into account 
intrinsic uncertainties arising from the ML technique 
(e.g., due to biases from the finite size of the event sam- 
ple or from an imperfect characterization of the CRB’s 
light curve), and are, most importantly, constructed to 
have proper coverage. Our final constraints on the LIV 
parameter and the LIV energy scale are produced using 
the calibrated CIs. 

As a final note, we would like to stress that the time 
shift T n (E n )i in Eq. (22) has been set to zero following 
Refs. [27, 29]. This implies that the time correction of any 
event entering the likelihood function is overestimated by 
a factor l/r] n , with 77^ = 1 — ( E n )i/E n £ [0.5, 1.0] for 
E £ [0.1, 30] GeV, n = 1 and, e.g., (E) t = 50 MeV. In 
principle, ignoring this time shift would thus produce an 
additional uncertainty r n — r n which is negative on aver- 
age. This would also slightly distort the likelihood func- 
tion since rj n varies with photon energy, possibly causing 
a reduction in sensitivity. In the large sample limit, one 
can show that the bias of the estimator takes the form 
b n ~ —T n (E n )i(E n )h/{E 2n )h, namely the fractional bias 
6 n /r n is negative and decreases with increasing hardness 
of the spectrum. In practice, it ranges from ^0.5% to 
~8% for spectral hardnesses similar to the ones of bursts 
we analyzed. In addition, due to the limited photon 
statistics available in our analysis and to the relatively 
small values of r n likely to be present in the data, the 
ratio b n /a[T n ] is also negligible (a few percent at most). 
One should, however, keep this effect in mind for future 
analyses of much brighter sources and/or in case of sig- 
nificant detections of LIV effects. 


C. Estimating the Systematic Uncertainty due to 
Intrinsic Spectral Evolution 


f{t) = C[t + T n (E n )i] / A; ~ C{t) / A,. (22) 

In practice, we fit the light curve C{t) with a function 
comprising up to three Gaussian functions (see for ex- 
ample Fig. 6). The fit is performed on events detected 
in a time interval somewhat wider than the default time 
intervals (defined in the beginning of Sec. IV) to allow 


So far we have concentrated on characterizing the sta- 
tistical uncertainties of our measurements. However, sys- 


5 These two values correspond to the (90%, 99%) CL quantile of 
a \i distribution. 
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tematic uncertainties can also be very important and 
should be taken into account, if possible. Here, we de- 
scribe how we model the dominant systematic uncer- 
tainty in our data, namely the intrinsic spectral evolu- 
tion observed in GRB prompt-emission light curves. A 
detailed discussion of this phenomenon is given in Sec. VI. 

For simplicity, in our introduction of LIV formalism 
(Sec. II), we implicitly ignored the presence of GRB- 
intrinsic effects (which can in general masquerade as a 
LIV-induced dispersion), and instead just used the quan- 
tity r n to describe the degree of dispersion in the data. 
However, r„ describes the total degree of dispersion and 
is, in general, the sum of the LIV-induced degree of dis- 
persion, described by a parameter tliv, and the GRB- 
intrinsic degree of dispersion, described by a parameter 
Tint, i.e., 

Tr, = Tint V Tliv . 

Our methods do not differentiate between the different 
sources of dispersion. Instead, they directly measure and 
constraint their sum r ra . We can either ignore any intrin- 
sic effects (i.e., assume Tj nt ~ 0) and proceed directly to 
constrain LIV using the obtained Cl on r n or we can first 
assume a model for ri nt , proceed to constrain tliv, and 
finally constrain LIV using the Cl on tliv- The second 
approach is more appropriate for constraining LIV, since 
its results are more robust with respect to the presence 
of GRB-intrinsic effects 6 . 

In principle, one could try to model r lnt using some 
knowledge of the physical processes generating the de- 
tected GRB emission or possibly using phenomenological 
models constructed from large sets of MeV / GeV observa- 
tions of GRBs. Unfortunately, because of the scarcity of 
GRB observations at LAT energies, neither approach has 
reached a mature enough stage to produce trustworthy 
and robust predictions of GRB spectral lags (at such en- 
ergies). Thus, any attempts to model Tj nt would, at this 
point, likely end up producing unreliable constraints on 
LIV. However, a more robust and conservative approach 
can be adopted, as follows. 

Since we do not have a model for Tj nt that reliably 
predicts GRB-intrinsic lags, we instead choose to model it 
in a way that produces the most reasonably conservative 
constraints on tliv- 

One of the main considerations behind modeling T; nt 
is the reasonable assumption that our measurements of 
r n are dominated by GRB intrinsic effects or in other 
words that our constraints on r„ also apply to r; nt . We 
start with the fact that we already have obtained a coarse 
measurement of the possible magnitude of T- lnt , provided 
by our constraints on r n . Specifically, we know that the 


Since the majority of previously published LIV constraints have 
not taken into account GRB-intrinsic effects, limits of the first 
approach are still useful for comparing experimental results 
across different studies. 


value of Ti n t (for a particular observation) is not likely 
larger than the width of the allowed range of r n , as de- 
scribed by its CIs 7 . Thus, we start with the working 
assumption that the width of the possible range of Tj nt is 
equal to the possible range of r n (as inferred by our CIs 
on it). 

Second, we assume that the T; nt has a zero value on 
average. This is a reasonable assumption given that we 
analyze in this study only cases where there is no clear 
detection of a spectral lag signal (i.e., t„ is consistent 
with zero within the uncertainty of its measured value). 
Moreover, this also avoids the need for introducing by 
hand a preferred sign for (Tj nt ). 

In principle, there are infinite choices for a particular 
shape of T- lnt given our constraint for its width and (zero) 
mean value. We choose the one that produces the least 
stringent (the most conservative) overall constraints on 
tliVj by modeling T; nt so that it reproduces the allowed 
range of possibilities of t„. For example, if our mea- 
surements imply that the data are compatible with (i.e., 
they cannot exclude) a positive r n , then we appropri- 
ately adjust Ti n t to match (explain) this possibility. This 
approach leads to confidence intervals on tliv that have 
the largest possible width. Other choices for modeling 
Tint can produce intervals more stringent either at their 
lower or their upper edge, but they cannot produce more 
stringent overall (i.e., when considering both their edges) 
constraints. The implementation of our model for Ti„ t , 
defined as TV in t(Rnt) with fi n t. being a random realiza- 
tion of T^t, depends on the particular method PV/SMM 
versus ML, and is described separately below. 

For constructing CIs on tliv with PV and SMM, we 
use a similar approach as in Eq. 11. However, instead of 
using as a pivotal quantity £ = f n — r n , we now use 

£' = r n — tliv 
= Tn T n T Tint 

— £ Rnt • ( 23 ) 

If we define the PDF of £' as Ps’ (£), where e' is a random 
realization of £' , and if q[i_cl)/2 an< ^ < t{i+CL)/2 are bs 
(1 — CL )/ 2 and (1 + CL )/ 2 quantiles, then starting from 
CL = Pr(q'^ 1 _ CL y 2 < £' < <^i +CL )/ 2 ), we derive a Cl on 
tliv of confidence level CL: 

CL = Pr[LL ' < tliv < UL') ( 24 ) 

= Pr{r n — q(i+cL)/2 < t liv < — Q(i-cl)/2)- 

Similarly to the Cl on r n which depends on the quan- 
tiles of Ps (approximated by / r ), the Cl on tliv depends 


7 The alternative case of a large Tint being approximately canceled 
by an oppositely large tliv seems extremely unlikely since it 
would require the improbable coincidence of LIV actually exist- 
ing, that the sign of the dispersion due to LIV being opposite of 
the sign of the dispersion due to intrinsic effects, and that the 
magnitudes of the two effects be comparable for each of the four 
GRBs (a “conspiracy of Nature”). 
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on the quantiles of Pgi. Assuming that the two compo- 
nents of £' (£ and Tint) are independent, Pg> is given by 
the convolution of their PDFs: 

/ OO 

/r(e' - f int )P Tint (f int )df int . (25) 

-CO 

Up to now, we have described a general way to pro- 
duce CIs on tliv, independent of the particular choice of 
P Ti nt ■ We mentioned above that we would like to choose 
a model for Tint such that it matches any part of the pa- 
rameter space for r n not excluded by the data. From the 
expression of the lower and upper limits for r n (Eq. 11) 
we observe that a large (say) positive tail in the f r distri- 
bution implies that our observations are compatible with 
(they cannot exclude) a symmetrically negative part of 
the parameter space of r„, and vice versa. Based on this 
observation, we choose to model P Tint (Tint) as f r (— Tint)- 
With this choice, Eq. 25 becomes: 

/ OO 

/r(e' - T int ) fr(- Tint) df int = Pac{ e'), 

-OO 

(26) 

where PacW) is defined as the autocorrelation func- 
tion of f r with argument e'. As an autocorrelation 
function, Pgi(t') is an even function with maximum at 
zero. Because of this symmetry, its (1 + CL)/ 2 and 
(1 - CL)/ 2 quantiles, q[ 1+CL y 2 and 9 (i-cl)/ 2> respec- 
tively, are equal. Thus, the confidence interval in Eq. 24 
is symmetric around the observed value f n . Finally, since 
(Tint) was chosen to be zero, in addition to f„ being our 
best estimate for t„, it is also our best estimate for tliv- 
A demonstration of the application of this method for 
GRB 090510, PairView, and n = 1 is shown in Sec. V, in 
Fig. 12. 

The confidence interval on tliv is wider than the one 
calculated on t„ by a degree that depends on the width 
and shape of the possible variations in T; nt (and thus of 
f r ). In the simple case of f r following a Gaussian dis- 
tribution, then the width would increase by a factor of 
\/2. In our case, the function f r does not always follow 
a Gaussian, hence the increase is not in general equal to 

V2. 

For the case of the ML method, we follow the same 
main idea (i.e., assume a P Tint following our observational 
uncertainty on t„ and produce confidence intervals on 
tliv) but apply it a different way. In this case, we run a 
second set of calibration simulations, in which the likeli- 
hood function is modified to include a not-necessarily- 
zero delay due to GRB-intrinsic effects. Specifically, 
Eq. 14 becomes: 

F'{E : £|tliv; Tint) = 

4 o E~ T e~ E / Et f(t - T LW E n - TintE n ) ■ (27) 

In each iteration of the simulation, we sample a differ- 
ent random value Tj nt from the assumed P Tint PDF and 
proceed normally to produce a distribution of lower and 
upper limits on tliv, the means of which will define our 


confidence interval on tliv- The P Tin t distribution is cho- 
sen in a similar way to the PV/SMM case using the distri- 
butions of r n produced during the first set of calibration 
simulations. The properties of the generated confidence 
intervals produced with this approach are the same as 
those constructed by the PV/SMM methods. 

We would like to add a point on the meaning of the 
distribution P Tint . In general, the properties of the emis- 
sion from a given GRB depend on two factors: the ini- 
tial properties describing the GRB’s generating system 
(e.g., mass, rotation speed, environment, redshift, etc.) 
and the randomness involved in the physical processes 
involved in producing the emission. We can imagine the 
Ti n t, quantity as a constant unknown parameter (a “true 
parameter”) that describes the range of possibilities for 
both factors mentioned above, thus P Tint can be consid- 
ered as its Bayesian prior. We can alternatively imagine 
the existence of some true parameter Ti ntj o = (Tint) (cho- 
sen to be zero) that depends solely on the progenitor 
properties, and that, during a GRB explosion, a random 
realization fi nt is produced depending on the Ti nt; o of that 
particular GRB system. In this case, we can imagine P Tint 
as a frequentist description of the range of possible Tint 
values occurring among an infinite number of GRBs, all 
initiated by the same initial conditions (i.e., having the 
same Ti n t,o). Based on the above, P rint can be consid- 
ered as a Bayesian prior or alternatively as a frequentist 
statement of the possible values of Tj nt across infinite rep- 
etitions of a GRB - the particular choice, however, does 
not matter. 

As a final note we should mention that our approach 
assumes that the experimental results allow the possibil- 
ity of T n being zero. With some additional assumptions, 
however, this approach can be generalized to include the 
case of a detection of a non-zero T n . For example, we 
could make the assumption that a detected non-zero to- 
tal dispersion is merely result of GRB-intrinsic effects, 
allow for (Ti n t) to take a non-zero value (with f n being 
the most conservative choice), and produce a final con- 
fidence interval on a residual tliv (that would still be 
consistent with a zero tliv)- 8 - It can be said that this 
method allows us to quantify the degree to which GRB- 
intrinsic effects reduce our ability to detect a residual 
LIV-induced dispersion. 

V. RESULTS 

The configuration of our methods is shown in Tab. II, 
in which we report the range (relative to the GBM trigger 


If a non-zero dispersion is detected, it would also be interesting to 
test the alternative possibility that this dispersion might indicate 
a non-zero value of tliv, rather than be fully attributed to ri n t 
as assumed in our method. Since most GRB properties vary 
weakly throughout the burst prompt emission, we may expect 
Tint t° also do so. In such a case, varying the time interval could 
change the measured value of Ti nt , while not affecting tliv- 
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time) of the analyzed data samples (common to all the 
methods), the value of SMM’s smoothing parameter p, 
the numbers of events used with PV and SMM Woo, and 
some quantities relevant to the ML method, namely the 
fitted index 7 of the observed spectrum S(E ), the number 
of events in the two parts of the data used for fitting 
the light-curve template N te mp iate and for calculating the 
likelihood 7Vfi t , and the energy separating these two parts 
of the data E cut . 

It is known that the spectra of LAT-detected GRBs 
typically comprise two spectral components: a Band (two 
smoothly connected power laws [49]) plus a power law 
function. These components do not necessarily have the 
exact same light curves and their spectra do not evolve 
in an identical fashion. As a result, an analysis of a data 
set consisting of events from both of these components 
might exhibit GRB-intrinsic spectral evolution that may 
be erroneously interpreted as LIV. This can be an im- 
portant systematic uncertainty, and is discussed further 
in Sec. VI. To reduce the influence of this effect, we per- 
formed the PV and SMM analyses on a data set starting 
from 100 MeV (instead of 30 MeV), a choice made a 
priori to reduce the contamination from the Band spec- 
tral component 910 . Because of the greater demand for 
statistics of the ML method, we did not apply such a 
minimum-energy cut for this method, and instead we 
used the events from E cut down to 30 MeV for the light- 
curve template construction. As a result, any differences 
in the temporal properties of the two spectral compo- 
nents might have affected the ML method more than the 
other two methods. However, the magnitudes of any such 
uncertainties are limited by the typically small contribu- 
tion of the Band component to the analyzed data and are 
likely smaller than the statistical errors. 

We produce constraints for two confidence levels: a 
90% two-sided (or equivalently 95% one-sided) CL and a 
99% two-sided (or equivalently 99.5% one-sided) CL. In 
the following, the “one-sided” or “two-sided” designations 
of the CLs may be omitted for brevity. 

An example plot used for choosing SMM’s p parame- 
ter, here for the case of GRB 090510 and n = 1, is shown 
in Fig. 4. For this case, we chose the value of p= 50, cor- 
responding to the minimum of the curve. The flatness 
of the curve around the minimum implies a weak depen- 
dence of the method’s sensitivity on p (in the vicinity of 
the minimum). 

Figure 5 demonstrates the application of the PV and 
SMM methods on GRB 090510 for n = 1. The top pan- 
els show how the best estimate of the LIV parameter is 


9 The spectrum of GRB 080916C comprises just one spectral com- 
ponent (Band). Thus, even though we did not need to reject the 
30—100 MeV events for that GRB, we still applied this cut for 
consistency between the four analyzed data sets. 

10 The particular value of 100 MeV is also the minimum energy typ- 
ically used in LAT science analyses, since the LAT reconstruction 
accuracy starts to deteriorate below this energy. 



FIG. 4. The median of the distribution of 99% CL upper 
limits (generated from simulated data sets inspired by the 
detected light curve) versus p. The error bars show the lo 
statistical uncertainty (arising from the finite number of sim- 
ulated data sets). This distribution is used for choosing the 
value of SMM’s p parameter for the GRB 090510 n = 1 ap- 
plication. 


measured, specifically from the location of the maximum 
of the KDE of the photon-pair lag distribution for PV 
(left column) and from the location of the maximum of 
the sharpness measure S for SMM (right column). The 
bottom panels show the distributions f r of the best LIV 
parameters in the randomized data sets, used for con- 
structing the CIs. Their asymmetry and features (in- 
versely) follow the shape of the analyzed light curves. 
The mean value of f r can be used as an estimate of the 
bias of T n ■ Except for GRB 090510, the magnitude of the 
bias is considerably smaller than the variance of f r (i.e. , 
up to ~ 10% of the variance); for GRB 090510, it in- 
creases up to 50% of the variance. The absolute value of 
the median of f r is for all cases smaller than ~ 10% of the 
variance. We correct f n for biases by subtracting from 
it the mean value of the f r distribution. The verifica- 
tion simulations of PV /SMM (described in Appendix A) 
show that the coverage of the produced CIs is approxi- 
mately proper even for asymmetric or non-zero-mean f r 
distributions, such as the ones shown. 

The light-curve template for GRB 090510 used by the 
ML method is shown in Fig. 6. Any statistical errors 
involved in the generation of the light-curve templates 
are properly included in the calibrated CIs of the ML 
method, as described in Appendix B. 

We show the spectral fit of the observed events from 
GRB 080916C in Fig. 7, which is used to calculate the 
spectral index 7 used by the ML method. The drop in 
the spectrum at low energies is caused by the sharp de- 
crease of the LAT effective area at those energies. In all 
cases, we choose E cut to be larger than the energy that 
this instrumental cutoff becomes important. This ensures 
that the spectral index 7 of the observed events is a good 
approximation of the index T of the incoming GRB flux 
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TABLE II. Configuration Details 


GRB Time Range (s) p N 100 7 -template N &t Scut (MeV) 

All Methods SMM PV & SMM Likelihood 



n = 1 

n — 2 

n = 1 

n = 2 

n = 1 

n = 2 

n ~ {1, 2} 

«={1,2} 

n — 1 

n = 2 

«={1,2} 

080916C 

3.53-7.89 

3.53-7.80 

50 

30 

59 

59 

2.2 

82 

59 

59 

100 

090510 

-0.01-3.11 

-0.01-4.82 

50 

70 

157 

168 

1.5 

148 

118 

125 

150 

090902B 

5.79-14.22 

5.79-14.21 

80 

80 

111 

111 

1.9 

57 

87 

87 

150 

090926A 

8.92-10.77 

9.3-10.76 

25 

30 

60 

58 

2.2 a 

53 

48 

47 

120 


a The spectral model for this GRB also includes an exponential cutoff with pre-set e-folding energy Tf =0.4 GeV in accordance with 
Ackermann et al. [48]. 



FIG. 5. Plots demonstrating the application of PV (left column) and SMM (right column) on GRB 090510 for n = 1. Top left: 
distribution of photon-pair lags (histogram), KDE of the distribution (thick curve), location of the KDE’s maximum used as 
r n by PV (vertical dashed line). Top right: sharpness measure S versus trial LIV parameter r„ (histogram), location of the 
curve’s maximum used as r n by SMM (vertical dashed line). Bottom row: distributions f r of the best estimates of the LIV 
parameter of the randomized data sets (histograms), 5% and 95% quantiles (dashed lines), 0.5% and 99.5% quantiles (dotted 
lines), and average value (central solid line). 


(within statistical errors). It allows us to considerably 
simplify the ML analysis by not having to deconvolve the 
instrument’s acceptance from the observed data or hav- 
ing to include the instrument’s response in the likelihood 
function. 

Finally, Fig. 8 demonstrates the application of the ML 
method, showing all the -2Aln(£) curves. We use the 
locations of the minima and the shapes of these curves 
to produce the best estimates and the (obtained directly 
from the data) CIs on r n , respectively. These curves 
are not exactly parabolic (and/or a transformation to a 


parabolic shape is not always possible). Therefore, any 
CIs produced based solely on their shape do not have 
an exactly proper coverage. The calibrated ML CIs (de- 
scribed in Appendix B) have by construction proper cov- 
erage, and are the ones used to constrain the quantities 
of the LIV models. 
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FIG. 6. Template light-curve fit for GRB 090510 used by 
the ML method. Actual data (histogram), fit template (thick 
curve), ±lcr error ranges of the fit (external thin curves). 



FIG. 7. Fit of the spectrum of detected events from 
GRB 080916C used to produce the index 7 of the S(E) func- 
tion used by the ML method. 


Constraints on the Total Degree of Dispersion, r n 

Table III reports our constraints on the total degree of 
dispersion, r n , and Fig. 9 shows our CIs on r„ plotted 
versus the distance n n . According to LIV models (i.e., 
Eq. 3), the magnitude of the observed dispersion due to 
LIV is proportional to K n . Thus, a positive correlation 
of r n and K n may imply a non-zero LIV effect. In our 
case and as can be seen from Fig. 9, no such correlation 
is evident. Additionally, all of our 99% CIs are compati- 
ble with a zero r n . Both features show that a LIV effect, 
if any, is dominated in this analysis by statistical and 
systematic (likely arising from GRB-intrinsic effects) un- 
certainties. Finally, we note that the results of the three 
methods (for the same GRB) are in good agreement to 




FIG. 8. -2Aln(£) curves for linear (top row) and quadratic 
(bottom row) dependence of the speed of light on energy. Left 
column: GRB 090510. Right column: GRBs 080916C (full 
line), 090902B (dotted line) and 090926A (dashed double- 
dotted line). 

each other (i.e., they have considerable overlap), evidence 
in support of the validity of each method. 

Table IV presents lower limits on Eqq calculated using 
our constraints on r„. The 95% lower limits are also 
plotted versus the redshift in Fig. 10. These limits do not 
take into account any GRB-intrinsic spectral evolution. 
Thus, while they are maximally constraining, they may 
not be as robust with regards to the presence of such 
intrinsic systematic uncertainties. 

Indeed, as we observed from, e.g., Fig. 10, some of our 
90% CL CIs are offset to a degree that their edges (i.e., 
limits) are very close to zero (e.g., GRB 090926A). For 
those CIs, the corresponding limits on Eqq are constrain- 
ing to a suspicious degree, given the considerably larger 
width of their CIs. It would be more acceptable if any 
very constraining limits were associated with correspond- 
ingly narrow CIs, contrary to what happens with some 
of the GRBs in our study. This feature required further 
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TABLE III. Our measurements on the LIV parameter r n describing the total degree of dispersion in the data. The limits are 
for a two-sided 99% CL. 


GRB Name PairView SMM Likelihood (from actual data) Likelihood (Calibrated) 3 

(Lower Limit, Best Value, Upper Limit) (sGeV^ 1 ) n = 1 


080916C 

-0.46 

0.69 

1.9 

-0.49 

0.79 

2.3 

-0.75 

0.1 

0.72 

-0.85 ■ 

0.77 

090510 (xlO 3 ) 

-73 

-14 

27 

-74 

-12 

30 

-25 

1 

6 

-9.8 

8.6 

090902B 

-0.36 

0.17 

0.53 

-0.25 

0.21 

0.62 

-0.25 

0.25 

0.55 

-0.63 ■ 

0.96 

090926A 

-0.45 

-0.17 

0.15 

-0.66 

-0.2 

0.23 

-0.45 

-0.18 

0.02 

-0.56 - 

0.18 





(Lower Limit, Best Value, Upper Limit) (sGeV 2 ) n 

= 2 


080916C 

-0.18 

0.45 

1.1 

-0.0031 

0.88 

2 

-0.9 

0.12 

1.1 

-0.83 

0.8 

090510 (xlO 3 ) 

-3.9 

-0.63 

0.88 

-4.1 

-0.68 

0.85 

-2.5 

-0.1 

0.3 

-0.32 

0.23 

090902B (xlO 3 ) 

-26 

17 

48 

-18 

24 

60 

-60 

10 

45 

-120 

110 

090926A 

-0.18 

-0.021 

0.13 

-0.12 

-0.06 

0.012 

-0.38 

-0.06 

0.11 

-0.44 - 

0.14 


a These are the ML CIs used for subsequently constraining LIV. 


TABLE IV. Lower Limits on Eqg for linear (n=l) and quadratic (n= 2) LIV for the subluminal (s± =+ 1) and superluminal 


(s±=-l) cases. 

The CL values 

are one-sided. These limits were produced using the total degree of dispersion 

in the data, r n . 

GRB Name 


PairView 

SMM 


Likelihood 3 




n=l, s±=+l (E 

pi units) 




95% 

99.5% 

95% 

99.5% 

95% 

99.5% 

080916C 

0.11 

0.081 

0.09 

0.067 

0.22 

0.2 

090510 

7.6 

1.3 

5.9 

1.2 

5.2 

4.2 

090902B 

0.17 

0.13 

0.15 

0.11 

0.12 

0.074 

090926A 

- 

0.55 

8 

0.35 

1.2 

0.45 




n=l, s±=-l (Ep; units) 




95% 

99.5% 

95% 

99.5% 

95% 

99.5% 

080916C 

18 

0.33 

5.4 

0.31 

0.2 

0.18 

090510 

0.56 

0.48 

0.57 

0.48 

11 

3.6 

090902B 

0.38 

0.2 

0.86 

0.28 

0.37 

0.11 

090926A 

0.24 

0.18 

0.2 

0.12 

0.17 

0.15 




n=2, s±=+l (10 10 

GeV un 

ts) 



95% 

99.5% 

95% 

99.5% 

95% 

99.5% 

080916C 

0.31 

0.28 

0.24 

0.21 

0.35 

0.33 

090510 

6.7 

3.3 

13 

3.3 

8.6 

6.4 

090902B 

0.8 

0.72 

0.73 

0.64 

0.64 

0.49 

090926A 

0.67 

0.48 

9.1 

1.6 

0.48 

0.47 




77=2, S±=-l (10 10 

GeV units) 



95% 

99.5% 

95% 

99.5% 

95% 

99.5% 

080916C 

- 

0.69 

- 

5.2 

0.34 

0.32 

090510 

1.9 

1.5 

1.9 

1.5 

9.4 

5.4 

090902B 

1.6 

0.97 

3.5 

1.2 

0.64 

0.46 

090926A 

0.51 

0.42 

0.51 

0.5 

0.31 

0.26 


a Calculated using the calibrated limits. 


scrutiny, hence, we examined our data and results in de- 
tail, and concluded that the CIs are offset likely because 
of GRB-intrinsic spectral evolution effects. 


For the case of GRB 090926A, the 90% CL CIs 
on T\ from our three methods, and the CIs on t 2 
from SMM for both CLs are either not consistent 
with zero or considerably offset towards negative val- 
ues (something which produces spuriously stringent up- 
per limits on r n ) . For example, the n = 1 CIs (not 


shown in the tables) are (-0.33, -0.17, -0.0010) s/GeV 11 
for PV, (-0.41, -0.20, 0.010) s/GeV for SMM, and 
(-0.25, -0.18, -0.13) s/GeV for ML (from data). As a re- 
sult, the 95% CL lower limits on Eqq i for the sublumi- 
nal case (s±=+l) are either suspiciously strong (SMM) 
or they could not be calculated at all (PV). The top 
left panel of Fig. 11 shows the E>100 MeV events from 
GRB 090926A processed by PV and SMM. As can be 
seen, the highest-energy photon in the data has an energy 


11 (lower limit, best estimate, upper limit) 
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FIG. 9. Our CIs on the total degree of dispersion in the data t „, obtained without taking into account any source-intrinsic 
effects, for linear (left panel) and quadratic (right panel) LIV. Each triplet of intervals corresponds to one GRB and shows, left 
to right, the results of PV, SMM, and ML (calibrated). The PV and ML points are drawn offset for visualization purposes. We 
present results for two CLs: a 90% (two-sided) CL denoted by the lines, and a 99% (two-sided) CL denoted by the external 
pairs of points. 


of ^3 GeV and is detected ^0.5 s before the main pulse. 
Our three methods predict that this event was most likely 
initially emitted in coincidence with the main pulse, and 
that it had been subsequently advanced by LIV to be de- 
tected before it. This case, shown in the top right panel 
of Fig. 11, implies a t\ ~-0.5s/3 GeV=-0.17 s/GeV, in 
accordance with the measured values. In the simulations 
performed for PV and SMM, such relatively small values 
were rare. Specifically, they occurred in a fraction of the 
iterations approximately equal to the ratio of the number 
of photons detected at least as early as the 3 GeV pho- 
ton (4) over the total number of photons (58 for n = 1), 
i.e. only 5-6%. This resulted in our 95% (one-sided) CL 
upper limits on n being negative or too small. 

The physical reason for these too negative CIs and 
Ti values may be GRB-intrinsic spectral evolution ef- 
fects, likely associated with the presence of spectral cut- 
off Ef ~0.4 GeV during the main bright pulse [48]. If 
this cutoff did not exist, more GeV photons might have 
been detected during this bright pulse, while if the cutoff 
also existed right before this pulse, the 3 GeV photon 
might have not been detected. Both cases would corre- 
spond to a fi closer to zero, and weaker, though, less 
spurious constraints. We conclude that our results from 
GRB 090926A are likely affected by a GRB-intrinsic spec- 
tral evolution, artificially strengthening (weakening) our 
limits on Eqq produced using r n for the subluminal (su- 
perluminal) case. 

Contrary to the case of GRB 090926A, for which the 
results hint towards negative t\ values, the results from 
GRB 080916C hint towards positive values. This either 
does not allow us to calculate lower limits on Eqq for the 
superluminal case (PV and SMM for n = 2; 95% CL) or 
produces spuriously constraining results (PV and SMM 
for n = 1 at 95% CL, and SMM n = 2 at 99.5% CL). 


A likely physical explanation for this positive lag is 
the progressive hardening of the prompt-emission spec- 
trum of GRB 080916C at LAT energies. According to 
broadband time- resolved spectroscopic studies [22], that 
spectrum can be adequately described by a Band func- 
tion, the high-energy component of which, /?, is initially 
very soft at a value of —2.63 ±0.12 during [0.004-3.58] s, 
hardens considerably to a value of —2.21 ± 0.03 during 
[3.58- 7.68] s, after which it stays constant (within statis- 
tics) to a value of —2.16 ± 0.03 up to at least 15.87 s. 
Based on this pattern, some soft-to-hard spectral evolu- 
tion is expected at least for the beginning of our analyzed 
intervals ([3.53-7.89] s for n = 1 and [3.53-7.80] s for 
n = 2). Similarly to the GRB 090926A case, we conclude 
that our GRB 080916C constraints on Eqq (produced 
using r n ) might also be affected by GRB spectral evolu- 
tion, artificially strengthening our superluminal-case lim- 
its and weakening our subluminal-case limits for PV and 
SMM. 


Finally, we notice that both of the calibrated ML lower 
limits on r ra for GRB 090510 are considerably more con- 
straining by about an order of magnitude than those from 
PV/SMM. We feel that this difference can be explained 
by the reduced sensitivity of the PV/SMM methods for 
constraining lower limits of the LIV parameter in the 
presence of long tails of the emission after the main peak, 
a feature of our chosen data set from GRB 090510. This 
effect was demonstrated in the one-to-one comparisons 
of the three methods described in Appendix C and illus- 
trated in the left panel of Fig. 20. Therefore, we attribute 
it to differences between the methods’ sensitivities. 
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FIG. 10. Our 95% (one-sided) CL lower limits on Eqg,i (left) and Uqg ,2 (right) LIV versus the redshift for s± = +1 (top; 
subluminal case) and s± = —1 (bottom; superluminal case). Similarly to Fig. 9, each triplet of markers corresponds to one 
GRB and shows limits calculated using the constraints on r n (i.e., without taking into account any source-intrinsic effects). 
The horizontal bars correspond to the averaged over the three methods lower limits on Eqg produced using the constraints 
on tliv (be., after accounting for GRB-intrinsic effects). On the left-hand plots we denote with the horizontal line the limit 
obtained by Fermi on GRB 090510 (DisCan; 95% limit obtained from paper’s Supplementary Information) [23]. On the top 
right plot we also denote the “high confidence” and “very high confidence” limits obtained by Fermi on GRB 090510 [23] and 
the 95% CL limit from H.E.S.S. study on PKS 2155-304 [29]. 


Constraints Using the LIV-Induced Degree of 
Dispersion, tliv 

The spuriously strong limits mentioned above imply 
that our sensitivity actually reaches the level of GRB- 
intrinsic effects. This motivated us to produce an addi- 
tional set of constraints, this time on tliv, taking into 
account intrinsic effects and according to the methodol- 
ogy in Sec. IV C. As an illustration of this method, Fig. 12 
shows the intermediate plots involved the calculation of 
the Cl on tliv for GRB 090510, PairView, and n = 1. 

For simplicity we do not report our CIs on tliv- In- 
stead, we just report the final limits on the LIV-model 
quantities, after averaging over the three methods. Ta- 
bles V and VI show our new 95% CL limits on Eqq 
and on the SME coefficients, respectively. Our lower lim- 


its on Eqg are also illustrated with the horizontal bars 
in Fig. 10, along with those produced without correct- 
ing for intrinsic effects (from T n ; shown with the mark- 
ers). As can be seen, the limits produced using tliv are 
considerably weaker than those produced using t„. The 
biggest difference is for the cases of GRBs 090926A and 
080916C, which had some spuriously strong limits that 
we attributed above to source-intrinsic effects. 


VI. SYSTEMATIC UNCERTAINTIES 

In this section, we discuss several systematic effects po- 
tentially influencing our results, namely those originating 
from the source and those having instrumental origins. 
Any dispersion induced by non-GRB standard physical 


18 


TABLE V. Our 95% CL lower limits on -Eqg, averaged over the three methods and calculated using the CIs on tliv (he., taking 
into account GRB-intrinsic effects). 


GRB Name 

s± = +1 

n = 1 ( Epi ) 

s± = -1 

n 

s± = +1 

= 2 (10 10 GeV) 

s± = -1 

080916C 

0.11 


0.32 

0.28 

0.56 

090510 

1.8 


3.2 

4.0 

3.0 

090902B 

0.11 


0.32 

0.58 

1.1 

090926A 

0.72 


0.15 

0.78 

0.41 


TABLE VI. Our 95% CL limits on the SME coefficients, averaged over the three methods and calculated using the CIs on tliv 
(i.e. , taking into account GRB-intrinsic effects). 


Model 

Source 

Quantity 

Lower Limit (10 20 GeV 2 ) 

Upper Limit (10“ 20 GeV” 2 ) 

Vacuum 

080916C 

y. 

0^(145% 120°)c^. m 

-8.7 

20 


090510 

y. 

oV/m(117 ,334 )c ( . J ^ m 

-0.31 

0.16 


090902B 

y. 

0 Y jm ( 63°,265°)c$. m 

-3.4 

5.2 


090926A 

y. 

0 Y jm ( 156°,353°)c$. m 

-11 

5.2 

Vacuum isotropic 

080916C 


(6) 

c moo 

-31 

70 


090510 


(6) 

c moo 

-1.1 

0.57 


090902B 


(6) 

c moo 

-12 

18 


090926A 


(6) 

C (J)00 

-37 

19 


processes is expected to be negligible compared to the 
dispersion produced by LIV [50] . 


Systematic Uncertainties from GRB-intrinsic Effects 

GRB-intrinsic effects that can cause systematic uncer- 
tainties in our results fall into two main categories: 

• the presence of multiple spectral components in the 
data not evolving with temporal coincidence, and 

• spectral evolution during the course of the burst or 
during each individual pulse. 

A full physical modeling of the emission processes occur- 
ring in the GRBs considered here is beyond the scope of 
this paper. Instead, we utilize published time-resolved 
spectral analyses to estimate the influence of any ob- 
served spectral evolution on our results. In the ini- 
tial Fermi papers on the GRBs analyzed in this study 
[22, 48, 51, 52], the prompt-emission spectra were fitted 
in relatively coarse time bins from keV to GeV energies 
with the combination of the empirical Band function with 
a high-energy power law. It was found that typically the 
Band component peaks at <MeV energies whereas the 
power-law component becomes dominant at LAT ener- 
gies (i.e., above ~100 MeV). 

In the case of GRB 080916C, the spectrum was well 
fitted by a Band function only, while the significance of 
the existence of an additional power-law component was 
found to be small. Some soft-to-hard spectral evolution 


could be present in the beginning of our analyzed in- 
tervals, as was discussed in the previous section. The 
broadband keV-GeV spectrum of the other three bursts 
is best represented by a combination of both spectral 
components: 

• in GRB 090510, the high-energy power law starts 
from the onset of the main emission in the LAT (at 
~0.5 s post-trigger) and dominates the Band com- 
ponent at energies above ~100MeV after ~0.7 s 
post-trigger. 

• In GRB 090902B, the high-energy power law is de- 
tected from the trigger time, and completely dom- 
inates the Band component in the LAT energy 
range. The spectral hardness of the emission in 
the LAT energy range is relatively constant during 
the time interval analyzed in this study. 

• In GRB 090926A, the high-energy power-law starts 
at the time of the bright pulse observed at ~10 s 
post-trigger and persists until ~22 s. Our analyzed 
time interval corresponds to the main bright pulse, 
during which the power law component dominates 
the emission in the LAT energy range, while ex- 
hibiting a high-energy spectral break with a cutoff 
energy Ef ~0.4 GeV. 

Since the two spectral components may be possibly 
originating from different physical regions of the burst 
and/or may be generated by physical processes evolving 
in different time scales, one might not necessarily expect 
them to be detected with exact temporal coincidence. 
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FIG. 11. SMM’s results for GRB 090926A and n = 1. 
Top left: photons detected in the default time interval, top 
right: SMM’s maximally-sharp version of these events. The 
curves/lines in the top row act as guides to the eye for the 
effects of a dispersion equal to the best measured value of fi=- 
0.17s/GeV (left figure) and zero (right figure). The bottom 
figures are of the same type as the figures in the right column 
of Fig. 5. 


This might lead to spurious signals originating from in- 
trinsic effects rather than LIV. There is only one case 
(GRB 090510) for which the LAT data during the ana- 
lyzed time intervals cannot be sufficiently approximated 
to contain a single spectral component, discussed in de- 
tail below. 

Using the spectral fits published in Ref. [51], we esti- 
mate that about half of the LAT-detected events from 
GRB 090510 below ~100-200MeV can be attributed to 
the Band component during the main episode observed 
around ~0.8 s post-trigger (comprising the bulk of the 
events in our analyzed time interval). This non-negligible 
fraction can potentially affect the ML method, which 
essentially compares the time profiles between the low 
(below ~100-150 MeV) and high energy emissions in 
the LAT. On the other hand, its effect on the PV and 
SMM methods is weaker because these methods ana- 
lyze a subset of the data produced with a higher-energy 
cut (E>100 MeV), for which only <15% of the events 
are estimated to be associated with the Band compo- 



FIG. 12. Demonstration of the generation of CIs on tliv for 
GRB 090510, PairView, and n = 1. The thin curve shows the 
normalized distribution f r used to approximate Ps{e) (the 
same distribution is also shown with different binning in the 
bottom left panel of Fig. 5), the thick curve shows the auto- 
correlation function P£/(e) calculated using Eq. 26, and the 
dashed lines show its 5% and 95% quantiles used for con- 
structing the 90% CL Cl on tliv- 


nent. Fortunately, evidence from cross checks performed 
in this work and from previously published results show 
that this effect has likely a negligible influence on our 
results. Specifically, a cross-correlation analysis between 
the time profiles of the keV-MeV emission (dominated 
by the Band component) and of the >100 MeV emission 
(dominated by the power-law component) of GRB 090510 
from 0.6 to 0.9 s [51] did not show any evidence of a time 
lag between the two spectral components. Furthermore, 
as shown in Appendix D, the PV and SMM CIs pro- 
duced using the data above 30 MeV are in good agree- 
ment with the results produced with the default cut of 
E > 100 MeV. We conclude that the inclusion of events 
related to the Band component for GRB 090510 did not 
cause any significant distortions in any of our analyses. 

Another potential source of systematic uncertainties 
is the spectral evolution detected with high significance 
in most LAT GRBs. One of its manifestations is the 
E>100 MeV emission detected by the LAT having a sys- 
tematically delayed onset with respect to the keV/MeV 
emission detected by the GBM [20]. Even though this 
spectral evolution can manifest as LIV, it happens so 
rapidly that typically only a very small fraction of the 
emission is detected during this transition. Furthermore, 
after the emission in the LAT is established, it usu- 
ally continues with a relatively stable degree of spectral 
hardness, at least according to the coarsely binned time- 
resolved spectral analyses mentioned above. 

For example, for the case of GRB 090510, cross- 
correlation analyses between the GBM-detected 
keV/MeV and LAT-detected E>100 MeV emissions re- 
vealed that the onset of the E>100 MeV emission trailed 
the onset of the keV/MeV emission by ~0.2-0.3 s [51]. 
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This offset implies the existence of a delay between 
the LAT data below and above 100 MeV, something 
that can potentially affect our results. However, the 
number of events detected during the onset of the LAT 
emission (~0.5~0.75 s) is negligible. Specifically, only 
~8% of the events above 30 MeV and within the default 
n = 1 interval were detected during the onset of the 
LAT emission. Furthermore, and as mentioned above, 
once the GRB 090510 emission is establishes in the LAT 
energy range, its spectral hardness remains relatively 
stable. We conclude that spectral evolution during the 
course of the emission of GRB 090510 affects only a very 
small fraction of the analyzed events. Thus, it is not 
expected to have a considerable influence on our results. 

GRB090926A is a peculiar case, as a strong spectral 
variability has been observed even after the onset of the 
high-energy emission in the LAT [48]. From the trigger 
time and until ~10 s, the high-energy power-law compo- 
nent is not detectable and the emission is well described 
by a single Band component. At ~10 s, a bright pulse 
appears (comprising the bulk of the events in our ana- 
lyzed time interval) , during which the power-law compo- 
nent becomes very bright dominating the emission and 
exhibiting a spectral cutoff at high energies. After the 
bright pulse, the two components become comparable in 
flux, while the cutoff of the power-law component ap- 
pears to be increasing in energy. Clearly, the results of 
a LIV analysis on an interval wide enough to include 
all these spectral-evolution effects would be strongly af- 
fected by them. By focusing only on a narrow snapshot 
of the GRB 090926 A’s emission (i.e. , the main bright 
pulse), during which the GRB spectrum is assumed not 
to vary too much 12 , we considerably reduced our expo- 
sure to such effects. 

At shorter time scales, the spectral hardness of GRB 
pulses is known to be correlated to their intensity and flu- 
ence at keV-MeV energies [55]. Due to the difficulty of 
measuring the GRB spectrum on a pulse-per-pulse basis 
with the limited photon statistics available to the LAT, 
there has been no evidence that this correlation extends 
to higher energies. However, the light curves of the GRBs 
analyzed in this study exhibit sharp peaks and fast vari- 
ability, thus the presence of any spectro-temporal cor- 
relations at high energies might, in principle, affect our 
results. This incomplete knowledge of GRB properties 
at high energies constitutes an intrinsic limitation of our 
model (e.g., it is unclear if the factorization in Eq. (14) 
holds at LAT energies at short time scales) and repre- 
sents a major source of systematic uncertainty in any 
GRB-based study of LIV. 


12 It should, however, be noted that even though an increase of 
the cutoff energy within the pulse could not be significantly de- 
tected due to the limited LAT statistics at GeV energies, an 
interpretation of this cutoff as arising from internal-opacity ef- 
fects does predict an associated evolution during the course of 
the spike [53, 54], 


Systematic Uncertainties from Instrumental Effects 

The probability of the LAT detecting an event of some 
energy depends on many factors, including the off-axis 
angle of the photon, with the probability being approx- 
imately inversely dependent on the off-axis angle. As a 
result, a constant-spectrum source observed at progres- 
sively larger (smaller) off-axis angles will correspond to 
a data set of a progressively harder (softer) average re- 
constructed energy. Such a data set may erroneously 
appear as containing a non-zero degree of spectral evolu- 
tion. Fortunately, this effect is negligible for our obser- 
vations since for the time scales under consideration the 
off-axis angles of the GRBs were almost constant. 

The energy-reconstruction accuracy of the LAT de- 
pends primarily on the true energies of the events. For 
the analyzed data set, about 90% of photons with energy 
greater than 1 GeV are predicted to have a reconstructed 
energy within ± ~20% of their true energy [40], which 
can be used to produce a rough estimate of the error 
on the produced limits on Eqq of up to 20% (90% CL). 
To verify this rough estimate we generated a collection 
of data sets derived from GRB 090510 by smearing the 
detected energies according to the energy dispersion func- 
tion of the instrument. For simplicity, during the produc- 
tion of the data sets we assumed that the detected energy 
was the true one. The 90% and 99% CL upper and lower 
limits varied by a factor of ~10% ( n = 1) and ~15% 
(n = 2) (1<t), in agreement with the rough estimate. 

The effective area of the LAT, corresponding to the 
P7_TRANSIENT_V6 selection used in this work, is typ- 
ically an increasing function of the energy up to ~ 
100 GeV. It starts from a zero value at few MeV and 
increases with increasing energy at a rate that is initially 
rapid but then gradually flattens above ~ 100 MeV. In 
the ML analysis, we have ignored the dependence of the 
effective area on the energy and approximated the spec- 
trum of incoming events with the spectrum of detected 
events (i.e., 7 ~ F). Because of this dependence, the 
spectrum of detected events appears slightly harder (less 
steep) than the spectrum of incoming events. This could 
affect the results of the ML analysis, depending on how 
sensitive it is on using an exactly correct spectral index. 
However, we have verified that the difference between the 
two spectral indices is always smaller than the statistical 
error of our measured spectral index, i.e., |y — T| < of. 
Thus, any systematic uncertainties by this approximation 
are dominated by the statistical uncertainty of determin- 
ing the true source spectrum. 

The effects from background contamination are ex- 
pected to be negligible, since the background rate for our 
data selection is very low, of the order of 0.1 (10 -3 ) Hz 
above 0.1 (1) GeV. 

The errors on the redshifts have a negligible effect on 
the lower limits on Eqq- Ala change in the redshift of 
GRB 080916C causes a relative change of about 10 -2 on 
the final limits. For the other GRBs in our sample, the 
relative change is also negligible, at the level of 10 -3 or 
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smaller. The errors on the cosmological parameters give 
an error of ~3%. 


VII. DISCUSSION/CONCLUSION 

We derive strong upper limits on the total degree of 
dispersion, r„, in the data of four LAT-detected GRBs. 
We use three statistical methods, one of which (PV) was 
developed as part of this study. The previously pub- 
lished most stringent limits on r n (at 95% CL; sub- 
luminal case) have been obtained for n = 1 by the 
Fermi GBM and LAT collaborations using GRB 090510 
{Eqg,i >3.5 .Epi; DisCan; Supplementary Information of 
Ref. [23]) 13 and for n = 2 by H.E.S.S. using the bright 
flares of PKS 2155-304 (£1qg,2 >6.4x 10 10 GeV; ML [29]). 
Our results from GRB 090510, namely -Eqg,i >7.6£'pi 
(PV) and -Eqg ,2 >1.3x 10 11 GeV (SMM), improve these 
constraints by a factor of ~2. 14 

In the above comparisons we do not consider other 
more constraining limits that were either produced in a 
very model-dependent manner or are of a moderate sta- 
tistical significance. Specifically, Chang et al. [56] tried to 
take into account intrinsic GRB time delays by estimat- 
ing them using the magnetic-jet GRB-emission model. 
However, our knowledge of GRB physics is not complete 
enough to be able to predict such intrinsic lags with suf- 
ficient certainty. Thus, even though such an approach 
proceeds in the right direction, it is highly sensitive to 
the particular choice and configuration of the employed 
model. Nemiroff et al. [25] took an innovative approach 
with which they zoomed in on the micro-structure of the 
burst’s emission above 1 GeV to produce very stringent 
constraints that were based, however, on observables of 
low statistical significance. 15 

To investigate why our GRB 090510 results are more 
constraining than the previous Fermi analysis of the 
same GRB, we applied the PV method to the same ex- 
act data used in the original Fermi publication. We used 
identical energy, time, and event selection cuts (as re- 
ported for the DisCan method), and obtained again more 


13 That work also reported lower limits on Eqq i as high as 10 Ep\. 
These limits, however, were not associated with a well quantified 
confidence level, but rather with a degree of confidence (“very 
high” to “medium”). Thus, they cannot be directly compared to 
the exact-CL limits produced in this work. 

14 At the 99% level, we improve on the Fermi limits Eqq i >1.2 Ep\ 
(DisCan) by a factor of ~ 4. 

15 They identified two pairs and one triplet of E>1 GeV photons in 
a 0.17 s interval of GRB 090510, with each photon being detected 
within ~1 ms of each other. The triplet, which contained the 
31 GeV photon, was used to place a stringent constraint. They 
calculated a probability of ~3 a for such a bunching of photons 
to arise by chance from a uniform emission in time. However, 
this significance is overestimated since it doesn’t account for the 
number of trials taken. Moreover, it does not reflect the con- 
fidence of their limit, since it strongly relies on associating the 
emission time of the 31 GeV photon with a tentative ms “pulse”. 


constraining results than the original Fermi publication 
by a factor of ~2-4 (depending on the CL). Addition- 
ally, we repeated our PV and SMM analyses using the 
configuration determined in this paper (i.e. , time inter- 
val, energy range, p) but using the P6_V3_TRANSIENT 
event selection of the previous Fermi work. The result- 
ing constraints were again of equal or higher strength (see 
Appendix D). These results show that the methods em- 
ployed in this work are more sensitive than the previous 
Fermi analyses. 

Our measurements are compatible with a zero degree 
of total dispersion in all the analyzed GRBs (at 99% CL). 
However, among these results, there are some spuriously 
strong limits on the total degree of dispersion, which we 
interpret as a product of GRB-intrinsic spectral-evolution 
effects. 

Using a maximally conservative set of assumptions 
to account for GRB-intrinsic effects, we constrain any 
residual dispersion in the data attributed solely to LIV, 
Tliv- The resulting CIs on tliv are less stringent than 
those on r n , albeit more robust with respect to the pres- 
ence of GRB-intrinsic effects, and thus, more appro- 
priate for constraining LIV. Our assumptions describe 
the worst-case scenario for GRB-intrinsic effects, and, as 
such, correspond to a maximum overall decrease in sen- 
sitivity. Our best constraints in the linear/subluminal 
case at 95% CL are Eqq,i >2 Ep\ for GRB 090510 and 
Eqg,i >0.1 Ep\ for the other three GRBs. We obtain re- 
sults of similar strength in the linear/super luminal case. 

As a final note we would like to mention that we con- 
sidered combining the results from the four GRBs to 
produce a single result that would be more constrain- 
ing and/or less affected by any GRB-intrinsic spectral- 
evolution effects (hoping that they might average out). 
However, we noticed that our GRB 090510 measurement 
is overall considerably more constraining than the other 
three cases. Thus, a combined result would not be very 
different from that of GRB 090510. Additionally, we do 
not expect that the intrinsic spectral-evolution effects for 
short GRBs (i.e., GRB 090510) are similar to those in 
long GRBs (other three cases). Thus, a combined analy- 
sis aimed at producing more robust results would have to 
be performed on short and long GRBs separately. Also, 
because our sample contains only a small number of long 
GRBs, we do not expect the average of their intrinsic 
effects to be an accurate representation of the typical 
long-GRB intrinsic evolution. Therefore, a combined re- 
sult obtained using the three long GRBs, would still be 
considerably less robust compared to each of the maxi- 
mally conservative CIs on tliv we produced here, e con- 
clude that there are no sufficient sensitivity or robustness 
benefits that a combined analysis of this limited data set 
can bring. 

There are many theoretical indications that Lorentz 
invariance may well break down at energies approaching 
the Planck scale. They come from the need to cut off 
the UV divergences in quantum field theory and black 
hole entropy calculations [57], from various quantum 
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gravity scenarios such as in loop quantum gravity [58], 
some string theory and M-theory scenarios, and non- 
commutative geometry models. There is one way to pre- 
scribe Lorentz invariance; there are many ways to violate 
Lorentz invariance. Kinematic tests of Lorentz invari- 
ance violation in QED depend on the possibility that 
the Lorentz violating terms can be different for electrons 
and photons [8] . It becomes even more complicated when 
hadronic interactions are considered. Many of these other 
tests, while quite sensitive, depend on the differences be- 
tween the individual maximum attainable velocities of 
various particle species [7]. In the context of effective 
field theories [35], birefringence tests have already pro- 
duced very strong constraints on LIV [10, 31]. Photo- 
hadronic interactions have also provided some powerful 
constraints [9]. 

One particular model inspired by string theory con- 
cepts presents the prospect that only photons would ex- 
hibit an energy-dependent velocity [32]. This model en- 
visions a universe filled with a gas of point-like D-branes 
that only interact with photons. It predicts that vacuum 
has an energy-dependent index of refraction that causes 
only a retardation. Since all photons are retarded, there 
is no associated vacuum birefringence effect, even though 
the degree of retardation has a first order dependence on 
the photon energy. The absence of an associated bire- 
fringence and the low-order ( n = 1) dependence of the 
predicted delay on the photon energy, render our results 
particularly unique for testing such a model 16 . Indeed, 
our constraints obtained using the total degree of disper- 
sion, t„, reiterate and support the previously-published 
results from Fermi [23], strongly disfavoring by almost 
an order of magnitude this model, and, in general, any 
class of models requiring Kqg.i 1$ Ep\- Our maximally- 
conservative set of constraints obtained using tliv also 
support the above statement. 

More GRB observations at high energies will allow 
us to improve GRB models and produce robust predic- 
tions on GRB-intrinsic delays (i.e., on P Tin t ), which will 
lead to more constraining CIs on tliv- Additionally, 
a larger collection of CIs on r n can be used for disen- 
tangling LIV-induced delays, which have a predicted de- 
pendence on the redshift, from the source-frame value 17 
of GRB-intrinsic delays, which can be assumed to not 


16 It has been argued that the D-brane model in Ref. [32] would 
suppress pair production interactions of ultrahigh energy (UHE) 
photons with cosmic microwave background photons, resulting 
in a flux of UHE photons in conflict with observations [59] . This 
would appear to be an independent argument against it. How- 
ever, in Ref. [60], it was argued that because electrons are not 
affected by the D-brane medium and because the pair produc- 
tion interaction involves an internal electron at the tree level, the 
resulting LIV effect in pair production is suppressed. Thus, this 
model is not ruled out by constraints on the UHE photon flux. 

17 The degree of intrinsic dispersion at the source is smaller than the 
observed degree of (intrinsic) dispersion at the Earth by a factor 
of (1 + z)™"*" 1 due to the relativistic expansion of the universe 
causing time dilation and redshift. 


have a redshift dependence or at least to have a differ- 
ent dependence than tliv ( see for example the approach 
in Refs. [16-18]). Future simultaneous observations of 
GRBs at MeV/GeV energies with Fermi-L AT and at 
GeV energies with HAWC [61] will have considerably in- 
creased statistics at GeV energies and a lever arm that 
extends to an even higher energy than this work, prop- 
erties that can provide uniquely constraining results. 
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Appendix A: PV and SMM Verification Tests 

We thoroughly tested the PV and SMM methods using 
a large number of simulated data sets to check for biases 
on the best estimates of the LIV parameter, to verify the 
proper coverage of the produced CIs, and to examine the 
robustness of the techniques (e.g., to find which proper- 
ties of the data could alter the validity and accuracy of 
the results). 

We performed the verification tests on a variety of col- 
lections of data sets, with each collection corresponding 
to a different light-curve and spectrum template, and to a 
different LIV parameter. The data sets of some collection 
represented the possible outcomes of the observation of 
the same exact source by a large number of identical de- 
tectors. By comparing the fraction of produced CIs that 
included the true value of the LIV parameter to their 
CL, we verified the coverage of these CIs. By repeating 
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this exercise on a diverse collection of data sets (pro- 
duced with, e.g., different statistics, number of pulses, 
light-curve asymmetry, pulse shape, spectral properties, 
degree of dispersion), we verified the robustness of the 
techniques. 

Our verification tests were performed on collections 
comprising ten thousand simulated data sets, with each 
of these sets being constructed in two steps: first its pho- 
ton energy-time pairs were randomly sampled from the 
light-curve and spectral template of the particular col- 
lection, and then a common degree of dispersion was ap- 
plied. 

We used two kinds of functional templates for the 
light curve. We started with simple synthetic templates 
composed of superpositions of Gaussian pulses of differ- 
ent widths, amplitudes, and means, and continued with 
templates inspired from the actually-observed GRB light 
curves. As an example, we show in Fig. 13 two of the 
light-curve templates used in our simulation. Both were 
inspired by actual detections, namely GRBs 090902B 
(top panel) and 090510 (bottom panel), representative 
cases of a long and short GRB, respectively. We obtained 
the light-curve templates using kernel density estimation 
(shown with the curve) of the actual light curves (his- 
tograms). 




FIG. 13. Two of the light-curve templates used in our ver- 
ification tests, one inspired by GRB 090902B (top) and one 
by GRB 090510 (bottom). The histograms are the GRB light 
curves (as actually detected), and the thick curves are our 
templates (produced by a KDE of the histograms). 

In all tested configurations, the energy spectrum of the 
non-dispersed data sets followed a power law, and ex- 
tended from 100 MeV to 40 GeV. For the GRB-based 
data sets we used an identical number of events as in the 
actual observations, and for the synthetic ones we simu- 
lated a range that was similar to that typically observed. 

We chose the maximum range of tested LIV parameters 
so that the simulated degree of dispersion did not distort 
the tested data too much. This way, we avoided the un- 
realistic possibility of having the highest-energy photons 
be disjoint and external from the bulk of the emission. 


To accomplish this, the magnitude of the tested LIV pa- 
rameter was not considerably larger than about the light- 
curve half-width divided by the highest simulated energy 
raised to the n power. 

We did not include any energy and temporal re- 
construction instrumental effects (i.e. , it was assumed 
that all photons were detected with the same energy- 
independent and constant-in-time efficiency) . A full sim- 
ulation including the LAT response to the GRB signal 
would also model any effects from a time-dependent effec- 
tive area and of any inaccuracies in the event-energy re- 
construction. The dependence of the results on both fac- 
tors is expected to be very small, as discussed in Sec. VI. 

As a demonstration of the verification process we 
present some of the diagnostic plots produced using the 
GRB 090510 light-curve template shown in the bottom 
panel of Fig. 13 and a zero LIV parameter. 
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FIG. 14. Top: 99% CL CIs produced by the application 
of PV on the GRB 090510-inspired simulated data set for 
n=0 s/GeV. The black dots denote the best estimates of the 
LIV parameter. Bottom: distributions of the lower (left) and 
upper limits (right) of these confidence intervals. The two 
external vertical lines denote the medians of the two distri- 
butions, and the middle vertical line denotes the mean of the 
best estimates. 

One of the first steps after a collection of data sets was 
constructed was to examine its distribution of associated 
confidence intervals. The top panel of Fig. 14 shows a 
stack of the confidence intervals produced by PV, and the 
bottom panel shows the distributions of lower and upper 
limits corresponding to these confidence intervals. The 
two external vertical lines in the latter figure denote the 
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distribution medians, and the middle vertical line shows 
the mean of the best estimates. By comparing the mean 
of the best estimates to the actual LIV parameter we 
checked for the presence of biases in the best estimates. 

Figure 15 shows two calibration plots produced by our 
simulations. These plots show the average best estimate 
and upper/lower limits on the LIV parameter for differ- 
ent injected values of r„. As can be seen, the methods 
properly measure the injected value with negligible bias. 
Furthermore, their sensitivity does not have a consider- 
able dependence on the injected degree of dispersion. 



FIG. 15. Calibration plots demonstrating some of the simu- 
lation results from the GRB 090902B- and 0900510-inspired 
data sets, top and bottom respectively. Each pair of intervals 
corresponds to a different value of the true (injected) LIV pa- 
rameter ti , and shows the results from PV (left interval) and 
SMM (right interval). The markers show the means of the 
best estimates, and the edges of the intervals correspond to 
the means of the upper and lower 99% CL limits on n. 


As mentioned in Sec. IV A, the distribution f r is used 
as an approximation of the PDF of the measurement er- 
ror of £, Pg. Since £ is a random variable (taking differ- 
ent values e across the simulated data sets), the quantity 
(7(e) = /p Pg(e)de is also a random variable. (7(e) be- 
haves similarly to a p- value, hence, C ~ (7(0, 1). We use 
the theoretical expectation of the uniformity of the PDF 
of (7, to verify whether the distribution f r (produced us- 
ing our randomization simulations described in Sec. IV A) 
is a good approximation of Pg , an approximation that is 
a cornerstone of our Cl-construction procedure. If this is 
indeed valid, then the empirical distribution, Pc emp , of 
the quantity C emp (ei ) = / Q e ‘ / rii (e)de, where and f r>i 
are the realizations of £ and f r in the i-th simulated data 
set, should also be distributed as a (7(0, 1). 

Figure 16 shows a normalized version of Pc Bmp pro- 
duced using the GRB 090510 inspired simulated data 
sets, PV, and t\ = 0 s/GeV. As can be seen, the empiri- 
cal distribution is indeed uniform, supporting the validity 
of our approximation f r ~ Pg. 

Pcemp is also used for verifying the coverage of the 
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FIG. 16. Normalized empirical distribution Pc Bmp (his- 
togram) obtained from the application of PV on the 
GRB 090510-inspired data set for tau i = 0 s/GeV super- 
imposed on its theoretical expectation (uniform distribution, 
shown with the horizontal line). The fact that the empiri- 
cal distribution follows a uniform distribution supports the 
validity of the assumptions behind our Cl construction. 


produced CIs, for any CL. 18 . Since the quantiles of the 
f r distribution are used for constructing our CIs, any 
erroneous distortions of f r (and equivalently any devia- 
tions of Pcemp f rom uniformity) will be associated with 
an improper coverage of the CIs. For example, if the 
CIs were erroneously under-covering, then Pc Bmp would 
acquire a V shape. On the other hand, if the CIs were 
erroneously wide (over covering), then Pc Bmp would ac- 
quire a A shape. By verifying the uniformity of Pc Bmp 
across its full range of values, we effectively tested the 
proper coverage of the CIs across the whole range of CLs 
(to the degree that the available statistics permitted). 

Using the verification tests mentioned above, we also 
found that 

• the sensitivities of both methods depend on the 
asymmetry on the light curve. Specifically, the 
longer the tail in the rising or falling side of the light 
curve is, the smaller the sensitivity of setting an up- 
per or lower limit on r„, respectively, becomes. The 
coverage, however, remains proper even for highly 
asymmetric light curves (e.g., like the one shown in 
the bottom panel of Fig. 13). 

• Miscoverage and bias can increase and sensitivity 
can decrease if the light curve includes separated 
bright pulses, due most likely to some form of in- 
terference between the individual pulses. This sys- 
tematic uncertainty becomes more prominent with 
the SMM method and when using large values of 
the p parameter. Our default data selection always 
includes a single bright pulse, so this problem does 
not affect our results. 

• Bias and miscoverage is larger for strongly spec- 
trally distorted light curves, i.e. , those produced 


18 We also performed the simple test of counting the fraction of 
CIs of a collection of data sets that included the true (injected) 
value of the LIV parameter to verify that the fraction was, as 
expected, equal to their CL, for two different values of CL: 90% 
and 99% 
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with a LIV parameter large enough that the 
highest-energy component is temporally disjoint 
from the bulk of the emission. The actual data did 
not appear to be spectrally distorted to the degree 
required for this systematic uncertainty to appear. 

• Tests performed on synthetic light curves compris- 
ing several pulses of a different spectral index (so as 
to simulate a GRB-intrinsic spectral evolution) re- 
vealed that this evolution is typically picked up by 
our methods as a non-zero LIV parameter, some- 
thing that reflects perhaps the most important ir- 
reducible uncertainty in our results. It is, however, 
fortunate that the additional GRB-intrinsic spec- 
tral evolution did not always dominate the simu- 
lation results, and that while there may be some 
non-zero bias, the miscoverage of the CIs was typ- 
ically not severe. 


Appendix B: Maximum Likelihood Method Tests 
and Calibrations 

Verification Tests 

We verified the ML method using Monte Carlo simu- 
lations, in a similar fashion to the PV/SMM methods, 
as described in the previous appendix. We performed 
tests on simple synthetic data sets as well as on data sets 
closely resembling the four GRBs in our sample. One of 
the main tests was the construction of calibration curves, 
in which we verified whether an injected LIV parameter 
was properly measured by the method with a reasonable 
degree of statistical accuracy. 

As an illustration of these tests, we show in Fig. 17 
a calibration plot demonstrating the simulation results 
from a GRB 0900510-inspired data set. The markers and 
the intervals show the average best estimates and 99% CL 
CIs on n, respectively. These averages were calculated 
across the different simulated realizations of the GRB 
emission. Our tests did not reveal any significant biases 
or other systematics. 



FIG. 17. A calibration plot obtained by the ML method on an 
GRB 090510-inspired data set. The markers show the means 
of the best estimates of n and the intervals correspond to the 
mean upper and lower 99% CL limits on n. 


Calibrated Confidence Intervals 

We construct calibrated CIs on r n by first generating 
several thousand simulated data sets having the same 
exact statistics as the actual data but with event energies 
and times randomly sampled from the fitted spectral and 
light-curve templates (e.g., such as from the templates 
shown in Figs. 7 and 6, respectively). We then apply 
the ML method to each one of them, using the same 
configuration as its application on the actual data, and 
calculate a Cl and a best estimate on r ra for each one 
of them. After all of the data sets have been processed, 
we calculate the average of the produced low and upper 
limits. Since we do not apply any spectral dispersion 
to the simulated data sets, we shift the mean low and 
upper limit values by the value of f n as measured from 
the actually detected data set, to finally produce a single 
calibrated CL 19 

The CIs are constructed using a pair of thresholds on 
— 2Aln(£) common to all the simulated data sets, and 
chosen to ensure the proper coverage of the produced CIs. 
Specifically, these thresholds are chosen so that exactly 
a fraction (1 — CL )/ 2 of the simulated lower limits and 
a fraction (1 + CL )/ 2 of the simulated upper limits are 
greater or smaller, respectively, than the value of r„ in 
the simulated data sets (equal by construction to zero). 

The calibration procedure includes the re-fitting of a 
light-curve template for each simulated sample. Thus, 
the produced CIs properly include the systematic uncer- 
tainties arising from the light-curve template generation 
procedure. On the other hand, for computational sim- 
plicity, we do not refit a spectral template, and instead 
use the one obtained from the actual data. Thus, the cal- 
ibrated CIs do not include uncertainties from the spectral 
fit. These are, however, negligible, since, as we have seen 
from toy Monte Carlo simulations and from the calcu- 
lations described in Sec. IV B, the final results depend 
weakly on the spectral index, contrary to their stronger 
dependence on the light-curve template. 

To illustrate the method, we show some intermediate 
results from its GRB 090510 application. Figure 18 shows 
the distributions of low and upper limits obtained from 
the simulated data sets. The mean values of these distri- 
butions are offset by t„ to produce our single calibrated 
upper and lower limits (i.e. , those shown in the last col- 
umn of Tab. III). From the mean of the simulated best 
estimates of the LIV parameter (see, e.g., Fig. 19) we esti- 
mated the bias of r n . In all cases, the bias was negligible 


19 In this last step and for simplicity, we make the assumption that 
the sensitivity of the method has a small dependence on r n , at 
least for the small possible values r n is expected to have (given 
past observations). Thus, we effectively assume that our sim- 
ulating a zero-LIV-parameter data set and then offsetting the 
mean upper and lower limits is equivalent to simulating a f n 
LIV-parameter data set and constructing a calibrated Cl directly 
from the mean lower and upper limits. 
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with respect to the root mean square of the simulated 
best estimates. 



FIG. 18. Distributions of the lower and upper 99% CL limits 
for n — 1 (left pad) and n = 2 (right pad) for GRB 090510. 
The vertical lines denote the means of the distributions, used 
for constructing the calibrated CIs. 
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FIG. 19. Distributions of the best LIV parameters obtained 
from each simulated data set for n = 1 (left) and n = 2 (right) 
for GRB 090510. Since there was no spectral dispersion ap- 
plied to the simulated data, these curves should peak near 
zero. 


For the ML method we used CIs calculated directly 
from the data (instead of from calibration simulations). 
However, we adjusted the two threshold values of - 
2Aln(£) used to produce its lower and upper limits, to 
ensure a proper coverage (evaluated across the simulated 
data set). 

The first and third panel of Fig. 20 show the obtained 
distributions of lower and upper limits, respectively. As 
can be seen, the sensitivities of the three methods are 
very similar. In the first panel, we also see that the ML 
method is slightly more sensitive when producing lower 
limits. We used this finding to explain in Sec. V why the 
ML method produced more constraining limits than the 
other two methods on n and GRB 090510. 

The histograms of the best estimates of the LIV param- 
eter (middle panel of Fig. 20) peak, as expected, near the 
true value of the LIV parameter, set equal to zero. The 
PV and SMM best-lag distributions peak at slightly more 
negative values than the approximately zero position of 
the ML method’s distribution. This can be attributed 
to the increased asymmetry of the PV and SMM distri- 
butions (skewness ~0.75) compared to the asymmetry 
of the likelihood distribution (skewness ^0.39), which 
moves the mode to lower values than the mean or the 
median. However, for considerations regarding the bias 
of the best estimates, the important fact is that both 
the median and the mean of these distributions are neg- 
ligible compared to their root mean square. Thus, the 
effect of any biases on the coverage of the produced CIs 
is expected to be negligible (as has been verified by the 
dedicated simulation tests). 

The 2D histograms in Fig. 21 provide a deeper view of 
how our methods compare. For the majority of the sim- 
ulated data sets, there is a close correspondence between 
their results. The PV and SMM results are the most 
similar, implying that these two methods probe the data 
in a similar fashion. The existence of differences between 
the methods’ results highlights their complementarity. 

Finally, we note that more than 99% of the exam- 
ined triplets of 90% CL CIs (one per simulated data set) 
are overlapping. This fraction is even larger (more than 
99.9%), if CIs of a higher CL (99%) are examined (not 
shown here). This large fraction of overlapping CIs shows 
that the troubling possibility of the three methods not 
allowing a common part of the parameter space is ex- 
tremely unlikely. 


Appendix C: Comparison of the Methods 

We compared the three methods by applying them to 
the same collection of simulated data sets to verify their 
validity and to help us explain any discrepancies observed 
in their application on the actual data. The simulated 
data of this test were produced using the GRB 090510- 
inspired light-curve template shown in Fig. 13 of Ap- 
pendix A and no extra applied dispersion (r„ was zero 
by construction). 


Appendix D: Analysis Cross-Checks 

We examined how the 99% CIs on r„ vary with respect 
to changes in the configuration of our methods and the 
data selection, to cross-check the validity and robustness 
of the results, and to gain insight on the behavior of our 
methods. Specifically: 

• we repeated the analysis excluding the highest- 
energy photon in the data, since it is expected to 
provide the most information about LIV dispersion. 
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FIG. 20. Histograms of the 95% (one-sided) CL lower limits (left panel), best estimates (middle panel), and upper limits (right 
panel) of the LIV parameter n produced by PV (black), SMM (red), and ML (blue) on simulated data sets. 
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FIG. 21. One to one comparisons of the 95% (one-sided) CL lower and upper limits (left and right columns respectively) and 
of the best estimates (middle column) on the LIV parameter n produced by our three methods on simulated data sets: ML vs 
SMM (top row), ML vs PV (middle row), SMM vs PV (bottom row). 


• We applied our methods on an extended time in- 
terval extending from the GRB trigger up to the 
time that the temporal variability has considerably 
subsided. The time intervals, selected with visual 
inspection, are 0-20 s for GRB 080916C, -0.01-10 s 
for GRB 090510, 0-60 s for GRB 090902B, and 0- 


40 s for GRB 090926A. The extended time intervals 
allow for maximal statistics, but at the same time 
potentially include a large degree of GRB-intrinsic 
spectral evolution that can, however, masquerade 
as LIV dispersion. 
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• we repeated the PV and SMM analyses using data 
produced using an earlier version of LAT’s event 
selection, P6_V3_TRANSIENT [21], also used by 
Fermi to constrain LIV [22, 23]. 

• Finally, we repeated the PV and SMM analy- 
ses starting from 30 MeV instead of their default 
100 MeV. While this change corresponds to in- 
creased statistics, it comes, however, with a larger 
contamination from the Band spectral component, 
which can increase the GRB-intrinsic systematic 
uncertainties. 

For the case of ML, we do not expect the calibrated 
CIs to vary in a considerably different way than the CIs 
obtained directly from the data, during the tests men- 
tioned above. Thus, for simplicity, we only present ML 
results obtained directly from the data. 

The test results are shown in Fig. 22. In all cases, 
the CIs produced by different methods (and for the same 
test) are in agreement with each other (i.e. , they have 
some overlap). Their widths and centers do change some- 
what across tests, something expected considering the 
different statistics and degrees of GRB-intrinsic spectral 
evolution in the different data sets. 

The removal of the highest- energy event, as expected, 
widened the produced CIs. The magnitude of the in- 
crease in their widths is a probe for the degree with which 
our methods draw information from the single highest- 
energy event and also for the systematic uncertainty as- 
sociated with the possibility of that event being back- 
ground. Because of the very low background contami- 
nation in our data sets, the highest-energy photons are 
typically securely associated to the GRB (see, e.g., the 
Supplementary Information of Ref. [23] regarding the as- 
sociation of the 31 GeV photon to GRB 090510). Thus, 
we do not consider the option of removing the highest- 
energy photon to increase the robustness of the results 
warranted. 

The changes brought by the use of the extended time 
interval did not correspond to a specific pattern. They 
were likely caused by the inclusion of emission of en- 
ergy considerably higher than that included in the de- 
fault time interval or the inclusion of significantly more 
GRB-intrinsic spectral evolution (likely in the case of 
GRB 090926 A). Perhaps the most significant change hap- 
pened with GRB 080916C and n = 2 on the PV and 
SMM results. For this case, the extended interval in- 
cluded at 13 GeV photon detected ~16.54 s post-trigger, 
which had an almost a decade in energy higher than that 
of the rest of the photons. As such, it dominated the 
PV/SMM estimation procedures with the edges of the 
confidence intervals on T 2 roughly corresponding to the 
time difference between its detection time and the edge of 
the analyzed interval divided by the square of its energy. 
The case of GRB 090926A is likely affected by both the 
inclusion of a very energetic event (~7 times higher en- 
ergy than the rest of the events) and the strong spectral 
evolution observed throughout this burst’s emission. We 


observe that the choice of time interval can significantly 
affect the final results, and conclude that an a priori and 
carefully chosen selection for the time interval, as in this 
work, is important for the validity of the results. 

Repeating the analysis with the P6_V3_TRANSIENT 
data set did not bring any considerable changes to the 
produced CIs, supporting the case that the improved lim- 
its produced in this work, when compared to past Fermi 
analyses of GRB 090510 [23], are a result of more sensi- 
tive analysis techniques rather than of a more constrain- 
ing data set. 

Finally, repeating the PV/SMM analyses starting from 
a lower minimum energy (30 MeV) did change the results 
significantly, implying that the systematic effects induced 
by the presence of two spectral components in the data 
are limited. 

We also repeated the ML analysis for GRBs 080916C 
and 090510 after performing some configurational 
changes affecting the light curve parametrization, such 
as using asymmetric (instead of symmetric) Gaussian 
pulses, a larger number of Gaussian pulses, different bin 
widths for the histogram used in producing the template 
(e.g., such as the one shown in Fig. 6), or different E cut 
values to split the data. The CIs varied up to a factor of 
~2 with respect to CIs obtained with the default config- 
uration. 
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FIG. 22. Comparison between the 99% CL CIs on r n ob- 
tained with the default configuration (horizontal dashed lines) 
to those obtained with modified configurations or data sets 
(vertical error bars). There are three groups of results in 
each figure, each corresponding to a different method: PV, 
SMM, ML (from left to right). The ML method CIs shown 
were produced directly from the data, instead using calibra- 
tion simulations. 
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